/*=========================================================================

  Program:   Visualization Toolkit
  Module:    vtkBoxClipDataSet.cxx
  Language:  C++
  Date:      $Date$
  Version:   $Revision$

  Copyright (c) 1993-2002 Ken Martin, Will Schroeder, Bill Lorensen 
  All rights reserved.
  See Copyright.txt or http://www.kitware.com/Copyright.htm for details.

  Copyright (C) 2003 Sandia Corporation
  Under the terms of Contract DE-AC04-94AL85000, there is a non-exclusive
  license for use of this work by or on behalf of the U.S. Government.
  Redistribution and use in source and binary forms, with or without
  modification, are permitted provided that this Notice and any statement
  of authorship are reproduced on all copies.

  Contact: Claudio Silva, csilva@cse.ogi.edu

     This software is distributed WITHOUT ANY WARRANTY; without even 
     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR 
     PURPOSE.  See the above copyright notice for more information.

=========================================================================*/
#include "vtkBoxClipDataSet.h"

#include "vtkCellArray.h"
#include "vtkCellData.h"
#include "vtkFloatArray.h"
#include "vtkGenericCell.h"
#include "vtkImageData.h"
#include "vtkIntArray.h"
#include "vtkMergePoints.h"
#include "vtkObjectFactory.h"
#include "vtkPointData.h"
#include "vtkUnsignedCharArray.h"
#include "vtkUnstructuredGrid.h"
#include "vtkIdList.h"

#include <math.h>
#include <stdio.h>

vtkCxxRevisionMacro(vtkBoxClipDataSet, "1.1.2.1");
vtkStandardNewMacro(vtkBoxClipDataSet);

//----------------------------------------------------------------------------
vtkBoxClipDataSet::vtkBoxClipDataSet()
{
  this->Locator = NULL;
  this->GenerateClipScalars = 0;

  this->MergeTolerance = 0.01;

  this->vtkSource::SetNthOutput(1,vtkUnstructuredGrid::New());
  this->Outputs[1]->Delete();
  this->InputScalarsSelection = NULL;
  this->Orientation = 1;
}

//----------------------------------------------------------------------------
vtkBoxClipDataSet::~vtkBoxClipDataSet()
{
  if ( this->Locator )
    {
    this->Locator->UnRegister(this);
    this->Locator = NULL;
    }
  this->SetInputScalarsSelection(NULL);
}

//----------------------------------------------------------------------------
// Do not say we have two outputs unless we are generating the clipped output. 
int vtkBoxClipDataSet::GetNumberOfOutputs()
{
  return 1;
}

//----------------------------------------------------------------------------
// Overload standard modified time function. If Clip functions is modified,
// then this object is modified as well.
unsigned long vtkBoxClipDataSet::GetMTime()
{
  unsigned long mTime=this->vtkDataSetToUnstructuredGridFilter::GetMTime();
  unsigned long time;

  if ( this->Locator != NULL )
    {
    time = this->Locator->GetMTime();
    mTime = ( time > mTime ? time : mTime );
    }

  return mTime;
}

//----------------------------------------------------------------------------
vtkUnstructuredGrid *vtkBoxClipDataSet::GetClippedOutput()
{
  if (this->NumberOfOutputs < 2)
    {
    return NULL;
    }
  
  return (vtkUnstructuredGrid *)(this->Outputs[1]);
}

//----------------------------------------------------------------------------
//
// Clip by box
//
void vtkBoxClipDataSet::Execute()
{
  vtkDataSet          *input = this->GetInput();
  vtkUnstructuredGrid *output = this->GetOutput();

  if (input == NULL)
    {
    return;
    }
  vtkIdType      i;
  vtkIdType      npts;
  vtkIdType     *pts;
  vtkIdType      estimatedSize;
  vtkIdType      newCellId;
  vtkIdType      numPts = input->GetNumberOfPoints();
  vtkIdType      numCells = input->GetNumberOfCells();
  vtkPointData  *inPD=input->GetPointData(), *outPD = output->GetPointData();
  vtkCellData   *inCD=input->GetCellData();
  vtkCellData   *outCD;
  vtkPoints     *newPoints;
  vtkPoints     *cellPts;
  vtkIdList     *cellIds;
  vtkIntArray   *locs;
  vtkDebugMacro(<< "Clip by Box");
  vtkUnsignedCharArray *types;

  int   j;
  int   cellType = 0;
  int   inputObjectType = input->GetDataObjectType();
  
  // if we have volumes
  if (inputObjectType == VTK_STRUCTURED_POINTS || 
      inputObjectType == VTK_IMAGE_DATA)
    {
    int dimension;
    int *dims = vtkImageData::SafeDownCast(input)->GetDimensions();
    for (dimension=3, i=0; i<3; i++)
      {
      if ( dims[i] <= 1 )
        {
        dimension--;
        }
      }
    }
  
  // Initialize self; create output objects
  //
  if ( numPts < 1 )
    {
    //vtkErrorMacro(<<"No data to clip");
    return;
    }
  
  // allocate the output and associated helper classes
  estimatedSize = numCells;
  estimatedSize = estimatedSize / 1024 * 1024; //multiple of 1024
  if (estimatedSize < 1024)
    {
    estimatedSize = 1024;
    }
  vtkCellArray *conn;
  conn = vtkCellArray::New();
  conn->Allocate(estimatedSize,estimatedSize/2);
  conn->InitTraversal();
  types = vtkUnsignedCharArray::New();
  types->Allocate(estimatedSize,estimatedSize/2);
  locs = vtkIntArray::New();
  locs->Allocate(estimatedSize,estimatedSize/2);
  newPoints = vtkPoints::New();
  newPoints->Allocate(numPts,numPts/2);
  
  // locator used to merge potentially duplicate points
  if ( this->Locator == NULL )
    {
    this->CreateDefaultLocator();
    }
  this->Locator->InitPointInsertion (newPoints, input->GetBounds());
  
  if ( !this->GenerateClipScalars && 
       !input->GetPointData()->GetScalars(this->InputScalarsSelection))
    {
    outPD->CopyScalarsOff();
    }
  else
    {
    outPD->CopyScalarsOn();
    }
  outPD->InterpolateAllocate(inPD,estimatedSize,estimatedSize/2);
  outCD = output->GetCellData();
  outCD->CopyAllocate(inCD,estimatedSize,estimatedSize/2);
  
  //Process all cells and clip each in turn

  vtkIdType       updateTime = numCells/20 + 1;  // update roughly every 5%
  vtkGenericCell *cell       = vtkGenericCell::New();

  int abort=0;
  int num;
  int numNew; 
  
  num       = 0;
  numNew    = 0;

  unsigned orientation = GetOrientation();   //Test if there is a transformation 

  //clock_t init_tmp = clock();
  for (vtkIdType cellId=0; cellId < numCells && !abort; cellId++)
    {
    if ( !(cellId % updateTime) )
      {
      this->UpdateProgress((float)cellId / numCells);
      abort = this->GetAbortExecute();
      }
    
    input->GetCell(cellId,cell);
    cellPts = cell->GetPoints();
    cellIds = cell->GetPointIds();
    npts = cellPts->GetNumberOfPoints();
                
    if (orientation)
      ClipHexahedron(newPoints,cell,this->Locator, conn,
         inPD, outPD, inCD, cellId, outCD);
    else
  
      ClipBox(newPoints,cell, this->Locator, conn,
          inPD, outPD, inCD, cellId, outCD);
    
    numNew = conn->GetNumberOfCells() - num;
    num    = conn->GetNumberOfCells();
    
    for (j=0; j < numNew; j++) 
      {
         locs->InsertNextValue(conn->GetTraversalLocation());
         conn->GetNextCell(npts,pts);
        
         //For each new cell added, got to set the type of the cell
         switch ( cell->GetCellDimension() )
           {
            case 0: //points are generated-------------------------------
              cellType = (npts > 1 ? VTK_POLY_VERTEX : VTK_VERTEX);
              break;
          
            case 1: //lines are generated----------------------------------
              cellType = (npts > 2 ? VTK_POLY_LINE : VTK_LINE);
              break;
          
            case 2: //polygons are generated------------------------------
              cellType = (npts == 3 ? VTK_TRIANGLE : 
                   (npts == 4 ? VTK_QUAD : VTK_POLYGON));
              break;
          
            case 3: //tetrahedra are generated------------------------------
              cellType = VTK_TETRA;
              break;
           } //switch
        
         newCellId = types->InsertNextValue(cellType);
         outCD->CopyData(inCD, cellId, newCellId);
      } //for each new cell
    } //for each cell
       
//  clock_t end_tmp = clock();
//  double  sec_tmp = ((double)end_tmp - init_tmp)/CLOCKS_PER_SEC;
//  fprintf(stdout," time : %f\n", sec_tmp);
  
  cell->Delete();
  
  output->SetPoints(newPoints);
  output->SetCells(types, locs, conn);

  conn->Delete();
  types->Delete();
  locs->Delete();
  
  newPoints->Delete();
  this->Locator->Initialize();//release any extra memory
  output->Squeeze();
}

//----------------------------------------------------------------------------
// Specify a spatial locator for merging points. By default, 
// an instance of vtkMergePoints is used.
void vtkBoxClipDataSet::SetLocator(vtkPointLocator *locator)
{
  if ( this->Locator == locator)
    {
    return;
    }
  
  if ( this->Locator )
    {
    this->Locator->UnRegister(this);
    this->Locator = NULL;
    }

  if ( locator )
    {
    locator->Register(this);
    }

  this->Locator = locator;
  this->Modified();
}

//----------------------------------------------------------------------------
void vtkBoxClipDataSet::CreateDefaultLocator()
{
  if ( this->Locator == NULL )
      {
    this->Locator = vtkMergePoints::New();
    this->Locator->Register(this);
    this->Locator->Delete();
      }
}

//----------------------------------------------------------------------------
// Set the box for clipping
// for each plane, specify the normal and one vertex on the plane. 
//
void vtkBoxClipDataSet::SetBoxClip(float *n0,float *o0,float *n1,float *o1,
                                   float *n2,float *o2,float *n3,float *o3,
                                   float *n4,float *o4,float *n5,float *o5)
{
  int i;

  for ( i=0;i<3;i++) 
    {
    n_pl[0][i] = n0[i];
    o_pl[0][i] = o0[i];
    }
  for ( i=0;i<3;i++) 
    {
    n_pl[1][i] = n1[i];
    o_pl[1][i] = o1[i];
    }
  for ( i=0;i<3;i++) 
    {
    n_pl[2][i] = n2[i];
    o_pl[2][i] = o2[i];
    }
  for ( i=0;i<3;i++) 
    {
    n_pl[3][i] = n3[i];
    o_pl[3][i] = o3[i];
    }
  for ( i=0;i<3;i++) 
    {
    n_pl[4][i] = n4[i];
    o_pl[4][i] = o4[i];
    }
  for ( i=0;i<3;i++) 
    {
    n_pl[5][i] = n5[i];
    o_pl[5][i] = o5[i];
    }
}

//----------------------------------------------------------------------------
// Specify the bounding box for clipping
void vtkBoxClipDataSet::SetBoxClip(float xmin,float xmax,
                                   float ymin,float ymax,
                                   float zmin,float zmax)
{
  SetOrientation(0);
  BoundBoxClip[0][0] = xmin;
  BoundBoxClip[0][1] = xmax;
  BoundBoxClip[1][0] = ymin;
  BoundBoxClip[1][1] = ymax;
  BoundBoxClip[2][0] = zmin;
  BoundBoxClip[2][1] = zmax;
}

//----------------------------------------------------------------------------
void vtkBoxClipDataSet::PrintSelf(ostream& os, vtkIndent indent)
{
  this->Superclass::PrintSelf(os,indent);

  os << indent << "Merge Tolerance: " << this->MergeTolerance << "\n";
  
  if ( this->Locator )
    {
    os << indent << "Locator: " << this->Locator << "\n";
    }
  else
    {
    os << indent << "Locator: (none)\n";
    }

  os << indent << "Generate Clip Scalars: " 
     << (this->GenerateClipScalars ? "On\n" : "Off\n");


  if (this->InputScalarsSelection)
    {
    os << indent << "InputScalarsSelection: " 
       << this->InputScalarsSelection << endl;
    }
}

//----------------------------------------------------------------------------
// TetraGrid: Subdivide cells in consistent tetrahedra.
// Case : Voxel(11) or Hexahedron(12).
//
// MinEdgF search the smallest vertex index in linear order of a face(4 vertices)
//
void vtkBoxClipDataSet::MinEdgeF(unsigned *id_v, vtkIdType *cellIds,unsigned *edgF)
{
    
  int i;
  unsigned int    id;
  int   ids;
  int min_f; 

  ids   = 0;
  id    = id_v[0];   //Face index     
  min_f = cellIds[id_v[0]];

  for(i=1;i<4;i++)
    {
    if(min_f > cellIds[id_v[i]])
      {
      min_f = cellIds[id_v[i]];
      id = id_v[i];            
      ids= i;
      }  
    }

  switch(ids)
    {
    case 0:
      if(id < id_v[2])
        {
        edgF[0] = id;
        edgF[1] = id_v[2];
        }
      else
        {
        edgF[0] = id_v[2];
        edgF[1] = id;
        }
      break;
    case 1:
      if(id < id_v[3])
        {
        edgF[0] = id;
        edgF[1] = id_v[3];
        }
      else
        {
        edgF[0] = id_v[3];
        edgF[1] = id;
        }
      break;
    case 2:
      if(id < id_v[0])
        {
        edgF[0] = id;
        edgF[1] = id_v[0];
        }
      else
        {
        edgF[0] = id_v[0];
        edgF[1] = id;
        }
      break;
  
    case 3:
      if(id < id_v[1])
        {
        edgF[0] = id;
        edgF[1] = id_v[1];
        }
      else
        {
        edgF[0] = id_v[1];
        edgF[1] = id;
        }
    break;
    }
}

//----------------------------------------------------------------------------
// TetraGrid: Subdivide cells in consistent tetrahedra.
//
// Case : Voxel or Hexahedron: 
//       if ( subdivide voxel  in 6 tetrahedra)
//           voxel : 2 wedges (*)
//       else 
//           voxel : 5 tetrahedra 
// 
// Case : Wedge (*)
// 
// ------------------------------------------------
//
//(*) WedgeToTetra: subdivide one wedge in  3 tetrahedra
//
//        wedge : 1 tetrahedron + 1 pyramid = 3 tetrahedra. 
//
void vtkBoxClipDataSet::WedgeToTetra(vtkIdType *wedgeId, vtkIdType *cellIds,vtkCellArray *newCellArray)
{
  int i;
  int id;
  vtkIdType xmin;
  vtkIdType newCellId;
  vtkIdType tab[4];
  vtkIdType tabpyram[5];

  static  vtkIdType vwedge[6][4]={{0,4,3,5},{1,4,3,5},{2,4,3,5},
                                          {3,0,1,2},{4,0,1,2},{5,0,1,2}};
        
  // the table 'vwedge' set 6 possibilities of the smallest index 
  //
  //             v5       
  //             /\       .
  //         v3 /..\ v4
  //           /   /
  //        v2/\  /
  //       v0/__\/v1
  //   if(v0 index ) is the smallest index:
  //   wedge is subdivided in:
  //         1 tetrahedron-> vwedge[0]: {v0,v4,v3,v5}
  //     and 1 pyramid    -> vert[0]  : {v1,v2,v5,v4,v0}
  //

  id = 0;     
  xmin = cellIds[wedgeId[0]];
  for(i=1;i<6;i++) 
    {
    if(xmin > cellIds[wedgeId[i]])
      {
      xmin = cellIds[wedgeId[i]];// the smallest global index
      id = i;                    // local index 
      }  
    }
  for (i =0;i<4;i++)
    {
    tab[i] = wedgeId[vwedge[id][i]];
    }
  newCellId = newCellArray->InsertNextCell(4,tab);
  
  // Pyramid :create 2 tetrahedra
  
  static vtkIdType vert[6][5]={{1,2,5,4,0},{2,0,3,5,1},{3,0,1,4,2},
               {1,2,5,4,3},{2,0,3,5,4},{3,0,1,4,5}};  
  for(i=0;i<5;i++)
    {
    tabpyram[i] = wedgeId[vert[id][i]];
    }
  PyramidToTetra(tabpyram,cellIds,newCellArray);
}

//----------------------------------------------------------------------------
// TetraGrid: Subdivide cells in consistent tetrahedra.
//
// PyramidToTetra :Subdivide the pyramid in consistent tetrahedra.
//        Pyramid : 2 tetrahedra.
//
void  vtkBoxClipDataSet::PyramidToTetra(vtkIdType *pyramId, vtkIdType *cellIds,vtkCellArray *newCellArray)
{
  vtkIdType xmin;
  vtkIdType newCellId;
  unsigned i,j,idpy;
  vtkIdType tab[4];

  // the table 'vpy' set  3  possibilities of the smallest index 
  // vertices{v0,v1,v2,v3,v4}. {v0,v1,v2,v3} is a square face of pyramid
  //
  //                v4      
  //                ^                
  //                                  
  //                       
  //           v3 _ _ __ _  v2     
  //           /         /     
  //        v0/_ _ _ _ _/v1         
  //                                  
  //   if(v0 index ) is the smallest index:
  //   the pyramid is subdivided in:
  //         2 tetrahedra-> vpy[0]: {v0,v1,v2,v4}
  //                        vpy[1]: {v0,v2,v3,v4}
  //    
  static vtkIdType vpy[8][4] ={{0,1,2,4},{0,2,3,4},{1,2,3,4},{1,3,0,4},
                                       {2,3,0,4},{2,0,1,4},{3,0,1,4},{3,1,2,4}}; 
                   
  xmin    = cellIds[pyramId[0]];
  idpy = 0;
  for(i=1;i<4;i++) 
    {
    if(xmin > cellIds[pyramId[i]])
      {
      xmin = cellIds[pyramId[i]];  // global index 
      idpy = i;                    // local index
      }
    }
  for(j = 0; j < 4 ; j++) 
    {
    tab[j] = pyramId[vpy[2*idpy][j]];
    }
  newCellId = newCellArray->InsertNextCell(4,tab);
  
  for(j = 0; j < 4 ; j++) 
    {
    tab[j] = pyramId[vpy[2*idpy+1][j]];
    }
  newCellId = newCellArray->InsertNextCell(4,tab);  
}


//----------------------------------------------------------------------------
//Tetra Grid : Subdivide cells in  consistent tetrahedra.
//             For each cell, search the smallest global index. 
//
//  Case Tetrahedron(10): Just insert this cell in the newCellArray
//
//  Case Voxel(11) or Hexahedron(12): 
//       - for each face: looking for the diagonal edge with the smallest index
//       - 2 possibilities: subdivide a cell in 5 or 6 tetrahedra
//
//         (I)Case 6 tetrahedra: 
//                  - 6 possibilities: subdivide de cell in 2 wedges: 
//                   (1) diagonal edges (v0,v5),(v2,v7) 
//                       vwedge[0]={0,5,4,2,7,6},
//                       vwedge[1]={0,1,5,2,3,7},
//                   (2) diagonal edges (v4,v7),(v0,v3)
//                       vwedge[2]={4,7,6,0,3,2},
//                       vwedge[3]={4,5,7,0,1,3},
//                   (3) diagonal edges (v0,v6),(v1,v7)
//        vwedge[4]= {1,7,5,0,6,4},
//                       vwedge[5]{1,3,7,0,2,6},
//                       subdivide de cell in 2 wedges:
//                   (4) diagonal edges (v1,v2),(v5,v6)
//                       vwedge[6]={4,5,6,0,1,2},
//                       vwedge[7]={6,5,7,2,1,3},
//                   (5) diagonal edges (v2,v4),(v3,v5)
//                       vwedge[8]={3,7,5,2,6,4},
//                       vwedge[9]={1,3,5,0,2,4},
//                   (6) diagonal edges (v1,v4),(v3,v6)
//                       vwedge[10]={0,1,4,2,3,6}
//                       vwedge[11]={1,5,4,3,7,6}
//         v6 _ _ __ _  v7     
//           /|        /|           VOXEL
//        v4/_|_ _ _ _/ |           opposite vertex of v0 is v7 and vice-versa
//          | |     v5| |           diagonal edges Edg_f[i]
//          |v2 _ _ _ |_|  v3        
//          |/        |/             
//        v0/_ _ _ _ _|/v1            
//
// 
//     (II)Case 5 tetrahedra: 
//       - search the smallest vertex vi
//       - verify if the opposite vertices of vi do not belong to any diagonal edges Edg_f
//       - 2 possibilites: create 5 tetraedra
//           - if vi is ( 0 or 3 or 5 or 6)  
//             vtetra[]={v0,v5,v3,v6},{v0,v4,v5,v6},{v0,v1,v3,v5},{v5,v3,v6,v7},{v0,v3,v2,v6}};      
//           - if vi is ( 1 or 2 or 4 or 7)                       
//             vtetra[]={v1,v2,v4,v7},{v0,v1,v2,v4},{v1,v4,v5,v7},{v1,v3,v2,v7},{v2,v6,v4,v7}};
//  Case Wedge (13): 
//       the table 'vwedge' set 6 possibilities of the smallest index 
//
//                 v5       
//               /\            .
//           v3 /..\ v4
//             /   /
//          v2/\  /
//         v0/__\/v1
//
//        if(v0 index ) is the smallest index:
//         wedge is subdivided in:
//              1 tetrahedron-> vwedge[0]: {v0,v4,v3,v5}
//               and 1 pyramid-> vert[0] : {v1,v2,v5,v4,v0}
//
//  Case Pyramid (14):
//       the table 'vpy' set  3  possibilities of the smallest index 
//       vertices{v0,v1,v2,v3,v4}. {v0,v1,v2,v3} is a square face of pyramid
//
//                v4  (opposite vertex of face with 4 vertices)  
//                ^                
//                                  
//                       
//          v3 _ _ __ _  v2     
//           /         /     
//        v0/_ _ _ _ _/v1         
//                                  
//        if(v0 index ) is the smallest index:
//        the pyramid is subdivided in:
//        2 tetrahedra-> vpyram[0]: {v0,v1,v2,v4}
//                       vpyram[1]: {v0,v2,v3,v4}
//  
//
//----------------------------------------------------------------------------
void vtkBoxClipDataSet::TetraGrid(vtkIdType typeobj, vtkIdType npts, vtkIdType *cellIds,vtkCellArray *newCellArray)
{
  vtkIdType tab[4];
  vtkIdType tabp[5];
  vtkIdType ptstetra = 4;
  vtkIdType newCellId;
  vtkIdType xmin;  
  unsigned i,j;
  unsigned id   =0;
  unsigned idpy =0;
        
  unsigned Edg_f[6][2]; //edge selected of each face
  unsigned idv[4];
    
  unsigned idopos;
  unsigned numbertetra;

  static  vtkIdType tetra[4]={0,1,2,3};
  static unsigned tabopos[8] = {6,7,4,5,2,3,0,1};

  switch(typeobj)
    {
    case VTK_TETRA:
      newCellId = newCellArray->InsertNextCell(ptstetra,tetra);
      break;  
    case VTK_VOXEL:
      // each face: search edge with smallest global index  
      // face 0 (0,1,5,4)
      idv[0]= 0;
      idv[1]= 1;
      idv[2]= 5;
      idv[3]= 4;
      MinEdgeF(idv,cellIds,Edg_f[0]);
      
      // face 1 (0,1,3,2)
      idv[0]= 0;
      idv[1]= 1;
      idv[2]= 3;
      idv[3]= 2;
      MinEdgeF(idv,cellIds,Edg_f[1]);
      // face 2 (0,2,6,4)
      idv[0]= 0;
      idv[1]= 2;
      idv[2]= 6;
      idv[3]= 4;
      MinEdgeF(idv,cellIds,Edg_f[2]);
      // face 3 (4,5,7,6)
      idv[0]= 4;
      idv[1]= 5;
      idv[2]= 7;
      idv[3]= 6;
      MinEdgeF(idv,cellIds,Edg_f[3]);
      // face 4 (2,3,7,6)
      idv[0]= 2;
      idv[1]= 3;
      idv[2]= 7;
      idv[3]= 6;
      MinEdgeF(idv,cellIds,Edg_f[4]);
      // face 5 (1,3,7,5)
      idv[0]= 1;
      idv[1]= 3;
      idv[2]= 7;
      idv[3]= 5;
      MinEdgeF(idv,cellIds,Edg_f[5]);
  
      //search the smallest global index of voxel 
      xmin = cellIds[0];
      id = 0;
      for(i=1;i<8;i++) 
        {
        if(xmin > cellIds[i])
          {
          xmin = cellIds[i];// the smallest global index
          id = i;           // local index 
          }  
        }
      //two cases: 
      idopos      = 7 - id;
      numbertetra = 5;
      for(i=0;i<6;i++) 
        {
        j=0;
        if (idopos == Edg_f[i][j]) 
          {
          numbertetra = 6;
          break;
          }
        j=1;
        if (idopos == Edg_f[i][j]) 
          {
          numbertetra = 6;
          break;
          }
        }
  
      if(numbertetra == 5) 
        {
        // case 1: create  5 tetraedra 
        if((id == 0)||(id == 3)||(id == 5)||(id == 6))
          {
          static  vtkIdType vtetra[5][4]={{0,5,3,6},{0,4,5,6},
                                      {0,1,3,5},{5,3,6,7},{0,3,2,6}}; 
          for(i=0; i<5; i++)
            {
            newCellId = newCellArray->InsertNextCell(4,vtetra[i]);
            }
          }
        else
          {
          static  vtkIdType vtetra[5][4]={{1,2,4,7},{0,1,2,4},
                                      {1,4,5,7}, {1,3,2,7},{2,6,4,7}};   
          
          for(i=0; i<5; i++)
            {
            newCellId = newCellArray->InsertNextCell(4,vtetra[i]);
            }
          }
        }
      else 
        {
        //case 2: create 2 wedges-> 6 tetrahedra
        static vtkIdType vwedge[12][6]={{0,5,4,2,7,6},{0,1,5,2,3,7},
                  {4,7,6,0,3,2},{4,5,7,0,1,3},
                  {1,7,5,0,6,4},{1,3,7,0,2,6},
                  {4,5,6,0,1,2},{6,5,7,2,1,3},
                  {3,7,5,2,6,4},{1,3,5,0,2,4},
                  {0,1,4,2,3,6},{1,5,4,3,7,6}}; 
        unsigned edgeId = 10*Edg_f[i][0]+ Edg_f[i][1];
        switch(edgeId)     
          {
          case 5:   // edge(v0,v5):10*0 + 5
          case 27:  // edge(v2,v7):10*2 + 7
            WedgeToTetra(vwedge[0],cellIds,newCellArray);
            WedgeToTetra(vwedge[1],cellIds,newCellArray);
            break;
          case 3:   // edge(v0,v2)
          case 47:  // edge(v4,v6)
            WedgeToTetra(vwedge[2],cellIds,newCellArray);
            WedgeToTetra(vwedge[3],cellIds,newCellArray);
            break;
          case 6:
          case 17:
            WedgeToTetra(vwedge[4],cellIds,newCellArray);
            WedgeToTetra(vwedge[5],cellIds,newCellArray);
            break;
          case 12:
          case 56:
            WedgeToTetra(vwedge[6],cellIds,newCellArray);
            WedgeToTetra(vwedge[7],cellIds,newCellArray);
            break;
          case 24:
          case 35:
            WedgeToTetra(vwedge[8],cellIds,newCellArray);
            WedgeToTetra(vwedge[9],cellIds,newCellArray);
            break;
          case 14:
          case 36:
            WedgeToTetra(vwedge[10],cellIds,newCellArray);
            WedgeToTetra(vwedge[11],cellIds,newCellArray);
            break;
          }
        }
        break;
    case VTK_HEXAHEDRON:
      // each face: search edge with smallest global index  
      // face 0 (0,1,5,4)
      idv[0]= 0;
      idv[1]= 1;
      idv[2]= 5;
      idv[3]= 4;
      MinEdgeF(idv,cellIds,Edg_f[0]);
      
      // face 1 (0,1,2,3)
      idv[0]= 0;
      idv[1]= 1;
      idv[2]= 2;
      idv[3]= 3;
      MinEdgeF(idv,cellIds,Edg_f[1]);
      // face 2 (0,3,7,4)
      idv[0]= 0;
      idv[1]= 3;
      idv[2]= 7;
      idv[3]= 4;
      MinEdgeF(idv,cellIds,Edg_f[2]);
      // face 3 (4,5,6,7)
      idv[0]= 4;
      idv[1]= 5;
      idv[2]= 6;
      idv[3]= 7;
      MinEdgeF(idv,cellIds,Edg_f[3]);
      // face 4 (3,2,6,7)
      idv[0]= 3;
      idv[1]= 2;
      idv[2]= 6;
      idv[3]= 7;
      MinEdgeF(idv,cellIds,Edg_f[4]);
      // face 5 (1,2,6,5)
      idv[0]= 1;
      idv[1]= 2;
      idv[2]= 6;
      idv[3]= 5;
      MinEdgeF(idv,cellIds,Edg_f[5]);
  
      //search the smallest global index of voxel 
      xmin = cellIds[0];
      id = 0;
      for(i=1;i<8;i++) 
        {
        if(xmin > cellIds[i])
          {
          xmin = cellIds[i];// the smallest global index
          id = i;           // local index 
          }  
        }
  
      //two cases: 
      idopos      = tabopos[id];
      numbertetra = 5;
      for(i=0;i<6;i++) 
        {
        j=0;
        if (idopos == Edg_f[i][j]) 
          {
          numbertetra = 6;
          break;
          }
        j=1;
        if (idopos == Edg_f[i][j]) 
          {
          numbertetra = 6;
          break;
          }
        }
  
      if(numbertetra == 5) 
        {
        // case 1: create  5 tetraedra 
        if((id == 0)||(id == 2)||(id == 5)||(id == 7))
          {
          static  vtkIdType vtetra[5][4]={{0,5,2,7},{0,4,5,7},
                                      {0,1,2,5},{5,2,7,6},{0,2,3,7}}; 
          for(i=0; i<5; i++)
            {
            newCellId = newCellArray->InsertNextCell(4,vtetra[i]);
            }
          }
        else
          {
          static  vtkIdType vtetra[5][4]={{1,3,4,6},{0,1,3,4},
                                           {1,4,5,6},{1,2,3,6},{3,7,4,6}};   
          
          for(i=0; i<5; i++)
            {
            newCellId = newCellArray->InsertNextCell(4,vtetra[i]);
            }
          }
        } 
      else 
        {
        //case 2: create 2 wedges-> 6 tetrahedra
        static vtkIdType vwedge[12][6]={{0,5,4,3,6,7},{0,1,5,3,2,6},
                  {4,6,7,0,2,3},{4,5,6,0,1,2},
                  {1,6,5,0,7,4},{1,2,6,0,3,7},
                  {4,5,7,0,1,3},{7,5,6,3,1,2},
                  {2,6,5,3,7,4},{1,2,5,0,3,4},
                  {0,1,4,3,2,7},{1,5,4,2,6,7}}; 
        unsigned edgeId = 10*Edg_f[i][0]+ Edg_f[i][1];

        switch(edgeId)
          {
          case 5:  // edge(v0,v5):10*0 + 5
          case 36: // edge(v3,v6):10*3 + 6
            WedgeToTetra(vwedge[0],cellIds,newCellArray);
            WedgeToTetra(vwedge[1],cellIds,newCellArray);
            break;
          case 2:  // edge(v0,v2)
          case 46: // edge(v4,v6)
            WedgeToTetra(vwedge[2],cellIds,newCellArray);
            WedgeToTetra(vwedge[3],cellIds,newCellArray);
            break;
          case 7:
          case 16:
            WedgeToTetra(vwedge[4],cellIds,newCellArray);
            WedgeToTetra(vwedge[5],cellIds,newCellArray);
            break;
          case 13:
          case 57:
            WedgeToTetra(vwedge[6],cellIds,newCellArray);
            WedgeToTetra(vwedge[7],cellIds,newCellArray);
            break;
          case 34:
          case 25:
            WedgeToTetra(vwedge[8],cellIds,newCellArray);
            WedgeToTetra(vwedge[9],cellIds,newCellArray);
            break;
          case 14:
          case 27:
            WedgeToTetra(vwedge[10],cellIds,newCellArray);
            WedgeToTetra(vwedge[11],cellIds,newCellArray);
            break;
          }
        }
      break;      
    case VTK_WEDGE:
      if(npts == 6) //create 3 tetrahedra
        {
        //first tetrahedron 
        static  vtkIdType vwedge[6][4]={{0,4,3,5},{1,4,3,5},{2,4,3,5},
                                       {3,0,1,2},{4,0,1,2},{5,0,1,2}};             
        xmin = cellIds[0];
        id = 0;
        for(i=1;i<6;i++) 
          {
          if(xmin > cellIds[i])
            {
            xmin = cellIds[i];// the smallest global index
            id = i;           // local index 
            }  
          }
        newCellId = newCellArray->InsertNextCell(4,vwedge[id]);
        
        //Pyramid :create 2 tetrahedra
           
        static vtkIdType vert[6][5]={{1,2,5,4,0},{2,0,3,5,1},{3,0,1,4,2},
                                   {1,2,5,4,3},{2,0,3,5,4},{3,0,1,4,5}};  
        static vtkIdType vpy[8][4] ={{0,1,2,4},{0,2,3,4},{1,2,3,4},{1,3,0,4},
                                   {2,3,0,4},{2,0,1,4},{3,0,1,4},{3,1,2,4}};                      
        xmin    = cellIds[vert[id][0]];
        tabp[0] = vert[id][0];
        idpy = 0;
        for(i=1;i<4;i++) 
          {
          tabp[i] = vert[id][i];
          if(xmin > cellIds[vert[id][i]])
            {
            xmin = cellIds[vert[id][i]]; // global index 
            idpy = i;                    // local index
            }
          }
        tabp[4] = cellIds[vert[id][4]];
        for(j = 0; j < 4 ; j++) 
          {
          tab[j] = tabp[vpy[2*idpy][j]];
          }
        newCellId = newCellArray->InsertNextCell(4,tab);
        
        for(j = 0; j < 4 ; j++) 
          {
          tab[j] = tabp[vpy[2*idpy+1][j]];
          }
        newCellId = newCellArray->InsertNextCell(4,tab);
        
        }
      else
        {
        fprintf(stdout," ERROR: this cell is not a wedge\n");
        }
      break;
  
    case 5: //VTK_PYRAMID nptos:5
                      //Create 2 tetrahedra
      if(npts== 5)
        {
        //note: the first element vpyram[][0] is the smallest index of pyramid
        static  vtkIdType vpyram[8][4]={{0,1,2,4},{0,2,3,4},{1,2,3,4},{1,3,0,4},
                  {2,3,0,4},{2,0,1,4},{3,0,1,4},{3,1,2,4}};
        xmin = cellIds[0];
        id   = 0;
        for(i=1;i<4;i++) 
          {
          if(xmin > cellIds[i])
            {
            xmin = cellIds[i]; // the smallest global index of square face 
            id = i;            // local index
            }          
          }
        newCellId = newCellArray->InsertNextCell(4,vpyram[2*id]);
        newCellId = newCellArray->InsertNextCell(4,vpyram[2*id+1]);
        }
      else
        {
        fprintf(stdout," ERROR: this cell is not a pyramid\n");
        exit(0);
        }
      break;
    }
}

//----------------------------------------------------------------------------
// The new cell  created in  intersection between tetrahedron and plane
// are tetrahedron or wedges or pyramides.
//
// CreateTetra is used to subdivide wedges and pyramides in tetrahedron.
//  
void vtkBoxClipDataSet::CreateTetra(vtkIdType npts,vtkIdType *cellIds,vtkCellArray *newCellArray)
{
  vtkIdType tabp[5];
  vtkIdType tab[3][4];
  
  unsigned i,j;
  unsigned id =0;
  unsigned idpy;

  vtkIdType xmin;
  vtkIdType newCellId;
        
  if (npts == 6) 
    { 
    //VTK_WEDGE: Create 3 tetrahedra
    //first tetrahedron 
    static  vtkIdType vwedge[6][4]={{0,4,3,5},{1,4,3,5},{2,4,3,5},
              {3,0,1,2},{4,0,1,2},{5,0,1,2}};             
    xmin = cellIds[0];
    id = 0;
    for(i=1;i<6;i++) 
      {
      if(xmin > cellIds[i])
        {
        xmin = cellIds[i];// the smallest global index
        id = i;           // local index 
        }  
      }     
    for(j = 0; j < 4 ; j++) 
      {
      tab[0][j] = cellIds[vwedge[id][j]];
      }  
    newCellId = newCellArray->InsertNextCell(4,tab[0]);
      
    //Pyramid: create 2 tetrahedra
  
    static vtkIdType vert[6][5]= {{1,2,5,4,0},{2,0,3,5,1},{3,0,1,4,2},
                  {1,2,5,4,3},{2,0,3,5,4},{3,0,1,4,5}};  
    static vtkIdType vpy[8][4] = {{0,1,2,4},{0,2,3,4},{1,2,3,4},{1,3,0,4},
                  {2,3,0,4},{2,0,1,4},{3,0,1,4},{3,1,2,4}};                      
    xmin    = cellIds[vert[id][0]];
    tabp[0] = vert[id][0];
    idpy = 0;
    for(i=1;i<4;i++) 
      {
      tabp[i] = vert[id][i];
      if(xmin > cellIds[vert[id][i]])
        {
        xmin = cellIds[vert[id][i]]; // global index 
        idpy = i;                    // local index
        }
      }
    tabp[4] = vert[id][4];
    for(j = 0; j < 4 ; j++) 
      {
      tab[1][j] = cellIds[tabp[vpy[2*idpy][j]]];
      }
    newCellId = newCellArray->InsertNextCell(4,tab[1]);
      
    for(j = 0; j < 4 ; j++) 
      {
      tab[2][j] = cellIds[tabp[vpy[2*idpy+1][j]]];    
      }
    newCellId = newCellArray->InsertNextCell(4,tab[2]);
      
    }
  else 
    {  
    //VTK_PYRAMID: Create 2 tetrahedra
    //The first element in each set is the smallest index of pyramid
    static  vtkIdType vpyram[8][4]={{0,1,2,4},{0,2,3,4},{1,2,3,4},{1,3,0,4},
                                      {2,3,0,4},{2,0,1,4},{3,0,1,4},{3,1,2,4}};
    xmin = cellIds[0];
    id   = 0;
    for(i=1;i<4;i++) 
      {
      if(xmin > cellIds[i])
        {
        xmin = cellIds[i]; // the smallest global index of face with 4 vertices 
        id = i;            // local index
        }          
      }
    for(j = 0; j < 4 ; j++) 
      {
      tab[0][j] = cellIds[vpyram[2*id][j]];
      }
    newCellId = newCellArray->InsertNextCell(4,tab[0]);
    
    for(j = 0; j < 4 ; j++) 
      {
      tab[1][j] = cellIds[vpyram[2*id+1][j]];
      }
    newCellId = newCellArray->InsertNextCell(4,tab[1]);
    } 
}

//----------------------------------------------------------------------------
// Clip each cell of an unstructured grid.
// 
//----------------------------------------------------------------------------   
//(1) How decide when the cell is NOT outside
//
//    Explaining with an example in 2D. 
//    Look at 9 region in the picture and the triangle represented there.
//    v0,v1,v2 are vertices of triangle T.
//
//              |         |        
//        1     |   2     |     3       
//     _ _ _ _ _ _ _ _ _  _ _ _ _ _ ymax
//              |         | v1        
//        4     |   5     |/\   6      
//              |         /  \             .
//     _ _ _ _ _ _ _ _ _ /| _ \_ _ _ymin
//              |     v2/_|_ _ \v0        
//        7     |   8     |     9  
//            xmin       xmax
//
// set test={1,1,1,1} (one test for each plane: test={xmin,xmax,ymin,ymax} )
//    for each vertex, if the test is true set 0 in the test table:
//    vO > xmin ?, v0 < xmax ?, v0 > ymin ?, v0 < ymax ?
//    In the example: test={0,1,1,0}
//    v1 > xmin ?, v1 < xmax ?, v1 > ymin ?, v1 < ymax ?
//    In the example: test={0,1,0,0}
//    v2 > xmin ?, v2 < xmax ?, v2 > ymin ?, v2 < ymax ?
//    In the Example: test={0,0,0,0} -> triangle is NOT outside.
//
//
//    In general, look at the possibilities of each region
//
//    (1,0,0,1)  | (0,0,0,1)  | (0,1,0,1)
//    - - - - - - - - - - - - - - - - - - 
//    (1,0,0,0)  | (0,0,0,0)  | (0,1,0,0)
//    - - - - - - - - - - - - - - - - - - 
//    (1,0,1,0)  | (0,0,1,0)  | (0,1,1,0)
//
//    In the case above we have:(v0,v1,v2)
//    (1,1,1,1)->(0,1,1,0)->(0,1,0,0)->(0,0,0,0)
//
//    If you have one vertex in region 5, this triangle is NOT outside
//    The triangle IS outside if (v0,v1,v2) are in region like 
//    (1,2,3): (1,0,0,1)->(0,0,0,1)->(0,0,0,1)
//    (1,4,7): (1,0,0,1)->(1,0,0,0)->(1,0,0,0)
//    (7,8,9): (1,0,1,0)->(0,0,1,0)->(0,0,1,0)
//    (9,6,3): (0,1,1,0)->(0,1,0,0)->(0,1,0,0) 
//
//    Note that if the triangle T is on the region 5, T is NOT outside
//    In fact, T is inside 
//
//    You can extend this idea to 3D.
//
//    Note: xmin = BoundBoxClip[0][0], xmax=  BoundBoxClip[0][1],... 
// 
//----------------------------------------------------------------------------  
// (2) Intersection between Tetrahedron and Plane:
//     Description: 
//         vertices of tetrahedron {v0,v1,v2,v3}
//         edge e1 : (v0,v1),   edge e2 : (v1,v2)
//         edge e3 : (v2,v0),   edge e4 : (v0,v3)
//         edge e5 : (v1,v3),   edge e6 : (v2,v3)
//         
//         interseting points p0, pi, ...
//
//         Note: The algorithm search intersection
//               points with these edge order.
//
//          v3                        v3
//       e6/|\ e5                     | 
//        / | \            e2         | 
//     v2/_ |_ \ v1    v2 - - -v1     |e4 
//       \  |  /                      |
//      e3\ | /e1                     |
//         \|/                        |
//          v0                        v0
//
//    (a) Intersecting  4 edges: see tab4[] 
//                    -> create: 2 wedges
//                    -> 3 cases
//                    ------------------------------------------
//        case 1246:   
//                   if (plane intersecting  edges {e1,e2,e4,e6})
//                   then  p0 belongs to e1, p1 belongs to e2,
//                         p2 belongs to e4, p3 belongs to e6.
//
//                                    v3
//                                     |            v3 
//                                     |            /
//                v1    v2 - *- -v1    |           * p3
//              /            p1        * p2       /
//             *p0                     |       v2/
//           /                         |
//          v0                        v0
//          (e1)          (e2)       (e4)        (e6) 
//
//                   So, there are two wedges:
//                     {p0,v1,p1,p2,v3,p3}  {p2,v0,p0,p3,v2,p1}
// 
//                   Note: if e1 and e2 are intersected by plane,
//                         and the plane intersects 4 edges,
//                         the edge e5 could not be intersected
//                         (skew lines, do not create a planar face)
//                         neither the edge e3((e1,e2,e3) is a face)
//
//                    ------------------------------------------
//        case 2345:   
//                   if (plane intersecting  edges {e2,e3,e4,e5})
//                   The two wedges are:
//                     {p3,v3,p2,p0,v2,p1}, {p1,v0,p2,p0,v1,p3}
//                         
//                    ------------------------------------------
//        case 1356:   
//                   if (plane intersecting  edges {e1,e3,e5,e6})
//                   The two wedges are:
//                     {p0,v0,p1,p2,v3,p3}, {p0,v1,p2,p1,v2,p3}
//
//                    -----------------------------------------
//
//    (b) Intersecting  3 edges: tab3[]
//                    ->create: 1 tetrahedron + 1 wedge
//                    -> 4 cases
//                    ------------------------------------------
//    case 134:
//                   if (plane intersecting  edges {e1,e3,e4})
//                   then  p0 belongs to e1, p1 belongs to e3,
//                         p2 belongs to e4.
//
//                           v3
//                           |             
//                           |            
//                     v2   *p2   v1      
//                       \   |   /  
//                      p1*  |  *p0     
//                          \|/
//                          v0
//                 
//
//                   So, there are:
//                       1 tetrahedron {v0,p0,p2,p1)
//                       and 1 wedge:    {p0,p2,p1,v1,v3,v2}: tab3[0] 
//                   
//                   Note:(a)if e1 and e3 are intersected by plane,
//                            and the plane intersects 3 edges,
//                            the edges e5 could not be intersected,
//                            if so, e6 be intersected too and the 
//                            plane intersect 4 edges.
//                        (b) tetrahedron vertices: 
//                            Use the first three indices of tab3:
//                              {v0,p(0),p(2),p(1)}                            
//                            
//                    ------------------------------------------  
//        case 125:
//                   if (plane intersecting  edges {e1,e2,e5})
//         There are :
//                      1 tetrahedron: {v1,p0,p2,p1}
//                      1 wedge      : {p0,p2,p1,v0,v3,v2},
//                    ------------------------------------------  
//    case 236:      
//                   if (plane intersecting  edges {e2,e3,e6})
//         There are :
//                      1 tetrahedron: {v2,p0,p1,p2}
//                      1 wedge      : {p0,p1,p2,v1,v0,v3},
//                    ------------------------------------------    
//        case 456:
//                  if (plane intersecting  edges {e4,e5,e6})
//         There are :
//                      1 tetrahedron: {v3,p0,p1,p2}
//                      1 wedge      : {p0,p1,p2,v0,v1,v2},
//                    ------------------------------------------
//
//    (c) Intersecting  2 edges and 1 vertex: tab2[] 
//                      -> create: 1 tetrahedron + 1 pyramid 
//      -> 12 cases 
//                    ------------------------------------------    
//        case 12:        v3        indexes:{0,0,1,2,3},
//                        *                  0 1 2 3 4      
//                    e6/ | \ e5                    
//                     /  |p1\         tetrahedron: indices(*,4,1,2)            
//                  v2/_ _|_*_\ v1     
//                    \   |   /                     
//                   e3\  |p0*                     
//                      \ | /e1                       
//                        v0                        
//
//                  if (plane intersecting  edges {e1,e2}
//                      and one vertex)
//           There are
//                      1 tetrahedron: {v1,v3,p0,p1}
//                      1 pyramid    : {v0,p0,p1,v2,v3},
//
//                  Note:  if e1 and e2 are intersected by plane,
//                            and the plane intersects 2 edges,
//                            then the vertex is v3.
//                    ------------------------------------------  
//      
//        case 13:                   indexes{2,1,0,1,3}
//                      Intersecting edges {e1,e3}
//                      Intersectinf vertes v3,
//           
//                      1 tetrahedron: {v0,v3,p1,p0}
//                      1 pyramid    : {v2,p1,p0,v1,v3},
//                    ------------------------------------------  
//        all cases see tab2[]
//                    ------------------------------------------
//    (d) Intersecting  1 edges and 2 vertices: tab1[] 
//                      -> create: 2 tetrahedra
//                      -> 6 cases
//
//        - edge e1, vertices {v2,v3}(e6)
//  
//                        v3        
//                        *              tetrahedra: {p0,v2,v3,v1},{p0,v3,v2,v0}
//                    e6/ | \ e5                     
//                     /  |  \                  .
//                 v2 *_ _|_ _\ v1     
//                    \   |   /                     
//                   e3\  |p0*                     
//                      \ | /e1                       
//                        v0                   
//         - other cases see tab1[]
//----------------------------------------------------------------------------  
//
void vtkBoxClipDataSet::ClipBox(vtkPoints *newPoints,
          vtkGenericCell *cell,
          vtkPointLocator *locator, 
          vtkCellArray *tets,
          vtkPointData *inPD, 
          vtkPointData *outPD,
          vtkCellData *inCD,
          vtkIdType cellId,
          vtkCellData *outCD)
{  
  vtkIdType     cellType   = cell->GetCellType();
  vtkIdList    *cellIds    = cell->GetPointIds();
  vtkCellArray *arraytetra = vtkCellArray::New();
  vtkPoints    *cellPts    = cell->GetPoints();
  vtkIdType     npts       = cellPts->GetNumberOfPoints(); 
  vtkIdType     cellptId[VTK_CELL_SIZE];
  vtkIdType     iid[4];
  vtkIdType    *v_id;
  vtkIdType    *verts, v1, v2;
  vtkIdType     ptId;
  vtkIdType     tab_id[6];
  vtkIdType     ptstetra = 4;

  vtkIdType i,j;
  unsigned allInside;
  unsigned ind[3];

  vtkIdType edges[6][2] = { {0,1}, {1,2}, {2,0},
                            {0,3}, {1,3}, {2,3} };  /* Edges Tetrahedron */
  float value,deltaScalar;
  float t, *p1, *p2;
  float v[3],x[3];
  float vf[4];
  float v_tetra[4][3];

  for (i=0; i<npts; i++)
    {
    cellptId[i] = cellIds->GetId(i);
    }

  // Convert all volume cells to tetrahedra       
  TetraGrid(cellType,npts,cellptId,arraytetra);   
  unsigned totalnewtetra = arraytetra->GetNumberOfCells();

  for (unsigned idtetranew = 0 ; idtetranew < totalnewtetra; idtetranew++) 
    {
    arraytetra->GetNextCell(ptstetra,v_id);

    // Test Outside: see(1)
    unsigned test[6] = {1,1,1,1,1,1};
    for (i=0; i<4; i++)
      {       
      cellPts->GetPoint(v_id[i],v);

      if (v[0] > BoundBoxClip[0][0])
        test[0] = 0;
      if (v[0] < BoundBoxClip[0][1])
        test[1] = 0;  
      if (v[1] > BoundBoxClip[1][0])
        test[2] = 0;
      if (v[1] < BoundBoxClip[1][1])
        test[3] = 0;
      if (v[2] > BoundBoxClip[2][0])
        test[4] = 0;
      if (v[2] < BoundBoxClip[2][1])
        test[5] = 0;
             
      }//for all points of the cell.

    if ((test[0] == 1)|| (test[1] == 1) ||
    (test[2] == 1)|| (test[3] == 1) ||
    (test[4] == 1)|| (test[5] == 1))
      {
        continue;                         // Tetrahedron is outside.
      }

    float  *pPtr = (float *)cellPts->GetVoidPointer(0);   
    for (allInside=1, i=0; i<4; i++)
      {
        ptId = cellIds->GetId(v_id[i]);
        cellPts->GetPoint(v_id[i],v);

        vf[0] = v[0];
        vf[1] = v[1];
        vf[2] = v[2];
        if (!(((vf[0] >= BoundBoxClip[0][0])&&(vf[0] <= BoundBoxClip[0][1]) &&
         (vf[1] >= BoundBoxClip[1][0])&&(vf[1] <= BoundBoxClip[1][1])&& 
         (vf[2] >= BoundBoxClip[2][0])&&(vf[2] <= BoundBoxClip[2][1]))))
          { 
          //outside,its type might change later (nearby intersection)
          allInside = 0;      
          }
    
        // Currently all points are injected because of the possibility 
        // of intersection point merging.
        if ( locator->InsertUniquePoint(v, iid[i]) )
          {
          outPD->CopyData(inPD,ptId, iid[i]);
          }
        
      }//for all points of the tetrahedron.

    if ( allInside )
      {   
        // Tetrahedron inside.
        vtkIdType newCellId = tets->InsertNextCell(4,iid);
        outCD->CopyData(inCD,cellId,newCellId);
        continue;
      }

    float  pc[3];
    float *pc1  , *pc2;
    float *pedg1,*pedg2;

      // tab4: tetrahedron intersection in 4 edges.see (2)

    unsigned tab4[6][6] = { {0,1,1,2,3,3},   // Tetrahedron Intersection Cases
              {2,0,0,3,2,1},
              {3,3,2,0,2,1},
              {1,0,2,0,1,3},
              {0,0,1,2,3,3},
              {0,1,2,1,2,3}};
    unsigned tab3[4][6] = { {0,2,1,1,3,2},
              {0,2,1,0,3,2},
              {0,1,2,1,0,3},
              {0,1,2,0,1,2}};
    unsigned tab2[12][5] = { {0,0,1,2,3},
               {2,1,0,1,3},
               {1,0,1,0,3},
               {1,0,1,2,0},
               {3,1,0,1,0},
               {2,0,1,3,0},
               {3,1,0,2,1},
               {2,1,0,0,1},
               {0,0,1,3,1},
               {1,0,1,3,2},
               {3,1,0,0,2},
               {0,0,1,1,2}}; 
    unsigned tab1[12][3] = { {2,3,1},
               {3,2,0},
               {3,0,1},
               {0,3,2},
               {1,3,0},
               {3,1,2},
               {2,1,0},
               {1,2,3},
               {0,2,3},
               {2,0,1},
               {0,1,3},
               {1,0,2}}; 
      
    vtkCellArray *cellarray = vtkCellArray::New();
    vtkIdType newCellId = cellarray->InsertNextCell(4,iid);

            // Test Cell intersection with each plane of box
    for (unsigned planes = 0; planes < 6; planes++) 
      {
      switch(planes)
        {
        case 0:
        case 1:
          ind[0] = 0;
          ind[1] = 1;
          ind[2] = 2;
          break;
        case 2:
        case 3:
          ind[0] = 1;
          ind[1] = 0;
          ind[2] = 2;
          break;
        case 4:
        case 5:
          ind[0] = 2;
          ind[1] = 0;
          ind[2] = 1;
          break;
        }
        
      value = BoundBoxClip[ind[0]][planes%2];   // The plane is always parallel to unitary cube. 
        
      unsigned totalnewcells = cellarray->GetNumberOfCells();
      vtkCellArray *newcellArray = vtkCellArray::New();
        
      for (unsigned idcellnew = 0 ; idcellnew < totalnewcells; idcellnew++) 
        {
        unsigned num_inter = 0;
        unsigned edges_inter = 0;
        unsigned i0,i1;
        float p_inter[4][3];   
        vtkIdType p_id[4];
        cellarray->GetNextCell(npts,v_id);
        
        newPoints->GetPoint(v_id[0],v_tetra[0]); //coord (x,y,z) 
        newPoints->GetPoint(v_id[1],v_tetra[1]); 
        newPoints->GetPoint(v_id[2],v_tetra[2]); 
        newPoints->GetPoint(v_id[3],v_tetra[3]); 

        for (int edgeNum=0; edgeNum < 6; edgeNum++)
          {
          verts = edges[edgeNum]; 
          
          p1 = v_tetra[verts[0]];
          p2 = v_tetra[verts[1]];
  
          if ( (p1[ind[0]] <= value && value <= p2[ind[0]]) || 
               (p2[ind[0]] <= value && value <= p1[ind[0]]) )
            {
            deltaScalar = p2[ind[0]] - p1[ind[0]];
      
            if (deltaScalar > 0)
              {
              pedg1 = p1;   pedg2 = p2;
              v1 = verts[0]; v2 = verts[1];
              }
            else
              {
              pedg1 = p2;   pedg2 = p1;
              v1 = verts[1]; v2 = verts[0];
              deltaScalar = -deltaScalar;
              }
      
            // linear interpolation
            t = ( deltaScalar == 0.0 ? 0.0 :
            (value - pedg1[ind[0]]) / deltaScalar );
      
            if ( t == 0.0 )
              {
              continue;
              }
            else if ( t == 1.0 )
              {
              continue;
              }
      
            pc1 = pPtr + 3*v1;
            pc2 = pPtr + 3*v2;
      
            for (j=0; j<3; j++)
              {
              x[j]  = pedg1[j]  + t*(pedg2[j] - pedg1[j]);
              pc[j] = pc1[j] + t*(pc2[j] - pc1[j]);
              p_inter[num_inter][j] = x[j];
              }
      
            // Incorporate point into output and interpolate edge data as necessary
            edges_inter = edges_inter * 10 + (edgeNum+1);
            if ( locator->InsertUniquePoint(x, p_id[num_inter]) )
              {
              outPD->InterpolateEdge(inPD, p_id[num_inter], cellIds->GetId(v1),
                   cellIds->GetId(v2), t);
              }
      
            num_inter++;        
            }//if edge intersects value
          }//for all edges

        if (num_inter == 0) 
          {           
          unsigned outside = 0;
          for(i=0; i<4; i++)
            {
            if (((v_tetra[i][ind[0]] < value) && ((planes % 2) == 0)) || 
                ((v_tetra[i][ind[0]] > value) && ((planes % 2) == 1))) 

              // If only one vertex is ouside, so the tetrahedron is outside 
              // because there is not intersection.
              {      
              outside = 1;
              break;
              }
            }
          if(outside == 0)
            {
             // else it is possible intersection if other plane
                 newCellId = newcellArray->InsertNextCell(4,v_id);       
            }
           continue;
          }
        switch(num_inter) 
          {
          case 4:                 // We have two wedges
            switch(edges_inter) 
              {
              case 1246:
                i0 = 0;
                break;
              case 2345:
                i0 = 2;
                break;
              case 1356:
                i0 = 4;
                break;
              default:
                printf(
                 "Intersection not found: Num_inter = %5d  Edges_inter = %8d\n",
                  num_inter,edges_inter);
              continue;
              }                                                         

            if (((v_tetra[3][ind[0]] < value) && ((planes % 2) == 0)) || 
                ((v_tetra[3][ind[0]] > value) && ((planes % 2) == 1))) 
              { 

              // The v_tetra[3] is outside box, so the first wedge is outside

              tab_id[0] = p_id[tab4[i0+1][0]];                            
              // ps: v_tetra[3] is always in first wedge (see tab)

              tab_id[1] = v_id[tab4[i0+1][1]];
              tab_id[2] = p_id[tab4[i0+1][2]];
              tab_id[3] = p_id[tab4[i0+1][3]];         
              tab_id[4] = v_id[tab4[i0+1][4]];
              tab_id[5] = p_id[tab4[i0+1][5]];
              CreateTetra(6,tab_id,newcellArray);
              }
            else 
              {
              tab_id[0] = p_id[tab4[i0][0]];
              tab_id[1] = v_id[tab4[i0][1]];
              tab_id[2] = p_id[tab4[i0][2]];
              tab_id[3] = p_id[tab4[i0][3]];
              tab_id[4] = v_id[tab4[i0][4]];
              tab_id[5] = p_id[tab4[i0][5]];
              CreateTetra(6,tab_id,newcellArray);
              }
            break;
          case 3:                   // We have one tetrahedron and one wedge
            switch(edges_inter) 
              {
              case 134:
                i0 = 0;
                break;
              case 125:
                i0 = 1;
                break;
              case 236:
                i0 = 2;
                break;
              case 456:
                i0 = 3;
                break;
              default:
                printf(
                "Intersection not found: Num_inter = %5d  Edges_inter = %8d\n",
                num_inter,edges_inter);
                continue;
              }                                                              

            if (((v_tetra[i0][ind[0]] < value) && ((planes % 2) == 0)) ||   
                ((v_tetra[i0][ind[0]] > value) && ((planes % 2) == 1))) 
              {

              // Isolate vertex is outside box, so the tetrahedron is outside
              tab_id[0] = p_id[tab3[i0][0]];
              tab_id[1] = p_id[tab3[i0][1]];
              tab_id[2] = p_id[tab3[i0][2]];
              tab_id[3] = v_id[tab3[i0][3]];
              tab_id[4] = v_id[tab3[i0][4]];
              tab_id[5] = v_id[tab3[i0][5]];
              CreateTetra(6,tab_id,newcellArray);
              }
            else 
              {
              tab_id[0] = v_id[i0];
              tab_id[1] = p_id[tab3[i0][0]];
              tab_id[2] = p_id[tab3[i0][1]];
              tab_id[3] = p_id[tab3[i0][2]];
              newCellId = newcellArray->InsertNextCell(4,tab_id);                
              }
            break;

          case 2:                    // We have one tetrahedron and one pyramid
            switch(edges_inter)      // i1 = vertex of the tetrahedron
              {
              case 12:
                i0 = 0; i1 = 1;  
                break;
              case 13:
                i0 = 1;  i1 = 0;
                break;
              case 23:
                i0 = 2;  i1 = 2; 
                break;
              case 25:
                i0 = 3;  i1 = 1;
                break;
              case 26:
                i0 = 4;  i1 = 2;
                break;
              case 56:
                i0 = 5;  i1 = 3;
                break;
              case 34:
                i0 = 6;  i1 = 0;
                break;
              case 46:
                i0 = 7;  i1 = 3;
                break;
              case 36:
                i0 = 8;  i1 = 2; 
                break;
              case 14:
                i0 = 9;  i1 = 0; 
                break;
              case 15:
                i0 = 10; i1 = 1;
                break;
              case 45:
                i0 = 11; i1 = 3;
                break;
              default:
                printf(
                    "Intersection not found: Num_inter = %5d  Edges_inter = %8d\n",
                  num_inter,edges_inter);
                continue;
              }                                                          

            if (((v_tetra[i1][ind[0]] < value) && ((planes % 2) == 0)) ||
                ((v_tetra[i1][ind[0]] > value) && ((planes % 2) == 1)))

              // Isolate vertex is outside box, so the tetrahedron is outside
              {
              tab_id[0] = v_id[tab2[i0][0]];
              tab_id[1] = p_id[tab2[i0][1]];
              tab_id[2] = p_id[tab2[i0][2]];
              tab_id[3] = v_id[tab2[i0][3]];
              tab_id[4] = v_id[tab2[i0][4]];
              CreateTetra(5,tab_id,newcellArray);
              }
            else 
              {
              tab_id[0] = v_id[i1];
              tab_id[1] = v_id[tab2[i0][4]];
              tab_id[2] = p_id[tab2[i0][1]];
              tab_id[3] = p_id[tab2[i0][2]];
              newCellId = newcellArray->InsertNextCell(4,tab_id);    
              }
            break;

          case 1:              // We have two tetrahedron.
            if ((edges_inter > 6) || (edges_inter < 1)) 
              {
              printf(
                "Intersection not found: Num_inter = %5d  Edges_inter = %8d\n",
                num_inter,edges_inter);
              continue;
              }                                                         
            if (((v_tetra[tab1[2*edges_inter-1][2]][ind[0]] < value) && 
                ((planes % 2) == 0)) ||   
                ((v_tetra[tab1[2*edges_inter-1][2]][ind[0]] > value) && 
                ((planes % 2) == 1))) 

              // Isolate vertex is outside box, so the tetrahedron is outside
              {   
              tab_id[0] = p_id[0];
              tab_id[1] = v_id[tab1[2*edges_inter-2][0]];
              tab_id[2] = v_id[tab1[2*edges_inter-2][1]];
              tab_id[3] = v_id[tab1[2*edges_inter-2][2]];
              newCellId = newcellArray->InsertNextCell(4,tab_id);
              }
            else 
              {
              tab_id[0] = p_id[0];
              tab_id[1] = v_id[tab1[2*edges_inter-1][0]];
              tab_id[2] = v_id[tab1[2*edges_inter-1][1]];
              tab_id[3] = v_id[tab1[2*edges_inter-1][2]];
              newCellId = newcellArray->InsertNextCell(4,tab_id);
              }
            break;
          }
        } // for all new cells

      cellarray->Delete();
      cellarray = newcellArray;
      } //for all planes

      unsigned totalnewcells = cellarray->GetNumberOfCells();

      for (unsigned idcellnew = 0 ; idcellnew < totalnewcells; idcellnew++) 
        {
        cellarray->GetNextCell(npts,v_id);
        newCellId = tets->InsertNextCell(npts,v_id);
        outCD->CopyData(inCD,cellId,newCellId);
        }
      cellarray->Delete();
    }
  arraytetra->Delete();
}

//----------------------------------------------------------------------------  
// ClipHexahedron: Box is like hexahedron.
// 
// The difference between ClipBox and ClipHexahedron is the outside test.
// The ClipHexahedron use plane equation to decide who is outside.
//
void vtkBoxClipDataSet::ClipHexahedron(vtkPoints *newPoints,
               vtkGenericCell *cell,
               vtkPointLocator *locator, 
               vtkCellArray *tets,
               vtkPointData *inPD, 
               vtkPointData *outPD,
               vtkCellData *inCD,
               vtkIdType cellId,
               vtkCellData *outCD) 
{  
  vtkIdType idcellnew;
  vtkIdType     cellType   = cell->GetCellType();
  vtkIdList    *cellIds    = cell->GetPointIds();
  vtkCellArray *arraytetra = vtkCellArray::New();
  vtkPoints    *cellPts    = cell->GetPoints();
  vtkIdType     npts       = cellPts->GetNumberOfPoints(); 
  vtkIdType     cellptId[VTK_CELL_SIZE];
  vtkIdType     iid[4];
  vtkIdType    *v_id;
  vtkIdType    *verts, v1, v2;
  vtkIdType     ptId;
  vtkIdType     tab_id[6];
  vtkIdType     ptstetra = 4;

  vtkIdType i,j;
  unsigned allInside;
  unsigned ind[3];

  vtkIdType edges[6][2] = { {0,1}, {1,2}, {2,0}, 
                            {0,3}, {1,3}, {2,3} };  /* Edges Tetrahedron */
  float deltaScalar;
  float t, *p1, *p2;
  float v[3],x[3];
  float p[6];
  float v_tetra[4][3];

  for (i=0; i<npts; i++)
    {
    cellptId[i] = cellIds->GetId(i);
    }
   
  TetraGrid(cellType,npts,cellptId,arraytetra);  // Convert all volume cells to tetrahedra

  unsigned totalnewtetra = arraytetra->GetNumberOfCells();

  for (unsigned idtetranew = 0 ; idtetranew < totalnewtetra; idtetranew++) 
    {
      arraytetra->GetNextCell(ptstetra,v_id);

      // Test Outside
      unsigned test[6] = {1,1,1,1,1,1};
      for (i=0; i<4; i++)
        {
        int k;
        cellPts->GetPoint(v_id[i],v);

        // Use plane equation 
        for(k=0;k<6;k++)
          {
          p[k] = n_pl[k][0]*
                (v[0] - o_pl[k][0])+ n_pl[k][1]*(v[1] - o_pl[k][1]) +  
                n_pl[k][2]*(v[2] - o_pl[k][2]);
          }

       
        for(k=0;k<3;k++)
          {
          if (p[2*k] < 0) test[2*k] = 0;
          if (p[2*k+1] < 0) test[2*k+1] = 0;
          }
             
        }//for all points of the cell.

      if ((test[0] == 1)|| (test[1] == 1) ||
          (test[2] == 1)|| (test[3] == 1) ||
          (test[4] == 1)|| (test[5] == 1))
        {
        continue;                         // Tetrahedron is outside.
        }
     
      float  *pPtr = (float *)cellPts->GetVoidPointer(0);   
      for (allInside=1, i=0; i<4; i++)
        {
        ptId = cellIds->GetId(v_id[i]);
        cellPts->GetPoint(v_id[i],v);
          
        for(unsigned k=0;k<6;k++)
            p[k] = n_pl[k][0]*(v[0] - o_pl[k][0])+ 
                   n_pl[k][1]*(v[1] - o_pl[k][1]) +  
                   n_pl[k][2]*(v[2] - o_pl[k][2]);

        if (!((p[0] <= 0) && (p[1] <= 0) &&
             (p[2] <= 0) && (p[3] <= 0) &&
             (p[4] <= 0) && (p[5] <= 0)))
          {
          allInside = 0;
          }
         
        // Currently all points are injected because of the possibility 
        // of intersection point merging.
        if ( locator->InsertUniquePoint(v, iid[i]) )
          {
          outPD->CopyData(inPD,ptId, iid[i]);
          }
        }//for all points of the tetrahedron.

      if ( allInside )
        {  
        vtkIdType newCellId = tets->InsertNextCell(4,iid);     // Tetrahedron inside.
        outCD->CopyData(inCD,cellId,newCellId);
        continue;
        }

      float  pc[3];
      float *pc1  , *pc2;
      float *pedg1,*pedg2;

      unsigned tab4[6][6] = { {0,1,1,2,3,3},   // Tetrahedron Intersection Cases
              {2,0,0,3,2,1},
              {3,3,2,0,2,1},
              {1,0,2,0,1,3},
              {0,0,1,2,3,3},
              {0,1,2,1,2,3}};
      unsigned tab3[4][6] = { {0,2,1,1,3,2},
              {0,2,1,0,3,2},
              {0,1,2,1,0,3},
              {0,1,2,0,1,2}};
      unsigned tab2[12][5] = { {0,0,1,2,3},
               {2,1,0,1,3},
               {1,0,1,0,3},
               {1,0,1,2,0},
               {3,1,0,1,0},
               {2,0,1,3,0},
               {3,1,0,2,1},
               {2,1,0,0,1},
               {0,0,1,3,1},
               {1,0,1,3,2},
               {3,1,0,0,2},
               {0,0,1,1,2}}; 
      unsigned tab1[12][3] = { {2,3,1},
               {3,2,0},
               {3,0,1},
               {0,3,2},
               {1,3,0},
               {3,1,2},
               {2,1,0},
               {1,2,3},
               {0,2,3},
               {2,0,1},
               {0,1,3},
               {1,0,2}}; 
      
      vtkCellArray *cellarray = vtkCellArray::New();
      vtkIdType newCellId = cellarray->InsertNextCell(4,iid);

      // Test Cell intersection with each plane of box
      for (unsigned planes = 0; planes < 6; planes++) 
        {
        switch(planes) 
          {
          case 0:
          case 1:
            ind[0] = 0;
            ind[1] = 1;
            ind[2] = 2;
            break;
          case 2:
          case 3:
            ind[0] = 1;
            ind[1] = 0;
            ind[2] = 2;
            break;
          case 4:
          case 5:
            ind[0] = 2;
            ind[1] = 0;
            ind[2] = 1;
            break;
          }
        
        vtkIdType totalnewcells = cellarray->GetNumberOfCells();
        vtkCellArray *newcellArray = vtkCellArray::New();
        
        for (idcellnew = 0 ; idcellnew < totalnewcells; idcellnew++) 
          {
          unsigned i0,i1;
          float p_inter[4][3];   
          unsigned num_inter = 0;
          unsigned edges_inter = 0;
          vtkIdType p_id[4];

          cellarray->GetNextCell(npts,v_id);
          
          newPoints->GetPoint(v_id[0],v_tetra[0]); //coord (x,y,z) 
          newPoints->GetPoint(v_id[1],v_tetra[1]); 
          newPoints->GetPoint(v_id[2],v_tetra[2]); 
          newPoints->GetPoint(v_id[3],v_tetra[3]); 

          p[0] = n_pl[planes][0]*(v_tetra[0][0] - o_pl[planes][0]) +
                 n_pl[planes][1]*(v_tetra[0][1] - o_pl[planes][1]) +  
                 n_pl[planes][2]*(v_tetra[0][2] - o_pl[planes][2]);  
          p[1] = n_pl[planes][0]*(v_tetra[1][0] - o_pl[planes][0]) + 
                 n_pl[planes][1]*(v_tetra[1][1] - o_pl[planes][1]) +  
                 n_pl[planes][2]*(v_tetra[1][2] - o_pl[planes][2]);
          p[2] = n_pl[planes][0]*(v_tetra[2][0] - o_pl[planes][0]) + 
                 n_pl[planes][1]*(v_tetra[2][1] - o_pl[planes][1]) +  
                 n_pl[planes][2]*(v_tetra[2][2] - o_pl[planes][2]);
          p[3] = n_pl[planes][0]*(v_tetra[3][0] - o_pl[planes][0]) + 
                 n_pl[planes][1]*(v_tetra[3][1] - o_pl[planes][1]) +  
                 n_pl[planes][2]*(v_tetra[3][2] - o_pl[planes][2]);
        
          for (int edgeNum=0; edgeNum < 6; edgeNum++)
            {
            verts = edges[edgeNum]; 
        
            p1 = v_tetra[verts[0]];
            p2 = v_tetra[verts[1]];
            float s1 = p[verts[0]];
            float s2 = p[verts[1]];    
            if ( (s1 * s2) <=0)
              {
              deltaScalar = s2 - s1;
            
              if (deltaScalar > 0)
                {
                pedg1 = p1;   pedg2 = p2;
                v1 = verts[0]; v2 = verts[1];
                }
              else
                {
                pedg1 = p2;   pedg2 = p1;
                v1 = verts[1]; v2 = verts[0];
                deltaScalar = -deltaScalar;
                t = s1; s1 = s2; s2 = t;
                }
        
              // linear interpolation
              t = ( deltaScalar == 0.0 ? 0.0 : ( - s1) / deltaScalar );
        
              if ( t == 0.0 )
                {
                continue;
                }
              else if ( t == 1.0 )
                {
                continue;
                }
        
              pc1 = pPtr + 3*v1;
              pc2 = pPtr + 3*v2;
        
              for (j=0; j<3; j++)
                {
                x[j]  = pedg1[j]  + t*(pedg2[j] - pedg1[j]);
                pc[j] = pc1[j] + t*(pc2[j] - pc1[j]);
                p_inter[num_inter][j] = x[j];
                }
        
              // Incorporate point into output and interpolate edge data as necessary
              edges_inter = edges_inter * 10 + (edgeNum+1);

              if ( locator->InsertUniquePoint(x, p_id[num_inter]) )
                {
                outPD->InterpolateEdge(inPD, p_id[num_inter], cellIds->GetId(v1),
                     cellIds->GetId(v2), t);
                }
        
              num_inter++;        
              }//if edge intersects value
            }//for all edges
          if (num_inter == 0) 
            {  
            unsigned outside = 0;
            for(i=0;i<4;i++)
              {
              if (p[i] > 0)
                {
                // If only one vertex is ouside, so the tetrahedron is outside 
                // because there is not intersection.
                // some vertex could be on plane, so you need to test all vertex

                outside = 1;  
                break; 
                }
              }
            if (outside == 0)
              {
              // else it is possible intersection if other plane

              newCellId = newcellArray->InsertNextCell(4,v_id);
              }
            continue;
            }
          switch(num_inter) 
            {
            case 4:                 // We have two wedges
              switch(edges_inter) 
                {
                case 1246:
                  i0 = 0;
                  break;
                case 2345:
                  i0 = 2;
                  break;
                case 1356:
                  i0 = 4;
                  break;
                default:
                  printf(
                    "Intersection not found: Num_inter = %5d  Edges_inter = %8d\n",
                    num_inter,edges_inter);
                  continue;
                }                                                    

              if (((p[3] > 0) && ((planes % 2) == 0)) ||
                  ((p[3] > 0) && ((planes % 2) == 1))) 
                {
                // The v_tetra[3] is outside box, so the first wedge is outside
                // ps: v_tetra[3] is always in first wedge (see tab)

                tab_id[0] = p_id[tab4[i0+1][0]];
                tab_id[1] = v_id[tab4[i0+1][1]];
                tab_id[2] = p_id[tab4[i0+1][2]];
                tab_id[3] = p_id[tab4[i0+1][3]];          
                tab_id[4] = v_id[tab4[i0+1][4]];
                tab_id[5] = p_id[tab4[i0+1][5]];
                CreateTetra(6,tab_id,newcellArray);
                }
              else 
                {
                tab_id[0] = p_id[tab4[i0][0]];
                tab_id[1] = v_id[tab4[i0][1]];
                tab_id[2] = p_id[tab4[i0][2]];
                tab_id[3] = p_id[tab4[i0][3]];
                tab_id[4] = v_id[tab4[i0][4]];
                tab_id[5] = p_id[tab4[i0][5]];
                CreateTetra(6,tab_id,newcellArray);
                }
              break;
            case 3:                   // We have one tetrahedron and one wedge
              switch(edges_inter) 
                {
                case 134:
                  i0 = 0;
                  break;
                case 125:
                  i0 = 1;
                  break;
                case 236:
                  i0 = 2;
                  break;
                case 456:
                  i0 = 3;
                  break;
                default:
                  printf(
                    "Intersection not found: Num_inter = %5d  Edges_inter = %8d\n",
                    num_inter,edges_inter);
                  continue;
                }                                                     

                if (((p[i0] > 0) && ((planes % 2) == 0)) ||   
                    ((p[i0] > 0) && ((planes % 2) == 1))) 
                  {
                  // Isolate vertex is outside box, so the tetrahedron is outside

                  tab_id[0] = p_id[tab3[i0][0]];
                  tab_id[1] = p_id[tab3[i0][1]];
                  tab_id[2] = p_id[tab3[i0][2]];
                  tab_id[3] = v_id[tab3[i0][3]];
                  tab_id[4] = v_id[tab3[i0][4]];
                  tab_id[5] = v_id[tab3[i0][5]];
                  CreateTetra(6,tab_id,newcellArray);
                  }
                else 
                  {
                  tab_id[0] = v_id[i0];
                  tab_id[1] = p_id[tab3[i0][0]];
                  tab_id[2] = p_id[tab3[i0][1]];
                  tab_id[3] = p_id[tab3[i0][2]];
                  newCellId = newcellArray->InsertNextCell(4,tab_id);                
                  }
                break;
            case 2:                    // We have one tetrahedron and one pyramid
              switch(edges_inter)      // i1 = vertex of the tetrahedron
                {
                case 12:
                  i0 = 0; i1 = 1;  
                  break;
                case 13:
                  i0 = 1;  i1 = 0;
                  break;
                case 23:
                  i0 = 2;  i1 = 2; 
                  break;
                case 25:
                  i0 = 3;  i1 = 1;
                  break;
                case 26:
                  i0 = 4;  i1 = 2;
                  break;
                case 56:
                  i0 = 5;  i1 = 3;
                  break;
                case 34:
                  i0 = 6;  i1 = 0;
                  break;
                case 46:
                  i0 = 7;  i1 = 3;
                  break;
                case 36:
                  i0 = 8;  i1 = 2; 
                  break;
                case 14:
                  i0 = 9;  i1 = 0; 
                  break;
                case 15:
                  i0 = 10; i1 = 1;
                  break;
                case 45:
                  i0 = 11; i1 = 3;
                  break;
                default:
                  printf(
                    "Intersection not found: Num_inter = %5d  Edges_inter = %8d\n",
                    num_inter,edges_inter);
                  continue;
                } 
              if (((p[i1] > 0) && ((planes % 2) == 0)) ||
                  ((p[i1] > 0) && ((planes % 2) == 1)))
                {
                // Isolate vertex is outside box, so the tetrahedron is outside

                tab_id[0] = v_id[tab2[i0][0]];
                tab_id[1] = p_id[tab2[i0][1]];
                tab_id[2] = p_id[tab2[i0][2]];
                tab_id[3] = v_id[tab2[i0][3]];
                tab_id[4] = v_id[tab2[i0][4]];
                CreateTetra(5,tab_id,newcellArray);
                }
              else 
                {
                tab_id[0] = v_id[i1];
                tab_id[1] = v_id[tab2[i0][4]];
                tab_id[2] = p_id[tab2[i0][1]];
                tab_id[3] = p_id[tab2[i0][2]];
                newCellId = newcellArray->InsertNextCell(4,tab_id);    
                }
              break;
            case 1:              // We have two tetrahedron.
              if ((edges_inter > 6) || (edges_inter < 1)) 
                {
                printf(
                  "Intersection not found: Num_inter = %5d  Edges_inter = %8d\n",
                  num_inter,edges_inter);
                continue;
                }                                        
              if (((p[tab1[2*edges_inter-1][2]] > 0) && ((planes % 2) == 0)) ||
                  ((p[tab1[2*edges_inter-1][2]] > 0) && ((planes % 2) == 1)))
                {
                // Isolate vertex is outside box, so the tetrahedron is outside
                tab_id[0] = p_id[0];
                tab_id[1] = v_id[tab1[2*edges_inter-2][0]];
                tab_id[2] = v_id[tab1[2*edges_inter-2][1]];
                tab_id[3] = v_id[tab1[2*edges_inter-2][2]];
                newCellId = newcellArray->InsertNextCell(4,tab_id);
                }
              else 
                {
                tab_id[0] = p_id[0];
                tab_id[1] = v_id[tab1[2*edges_inter-1][0]];
                tab_id[2] = v_id[tab1[2*edges_inter-1][1]];
                tab_id[3] = v_id[tab1[2*edges_inter-1][2]];
                newCellId = newcellArray->InsertNextCell(4,tab_id);
                }
              break;
            }
          } // for all new cells
        cellarray->Delete();
        cellarray = newcellArray;
        } //for all planes

      vtkIdType totalnewcells = cellarray->GetNumberOfCells();

      for (idcellnew = 0 ; idcellnew < totalnewcells; idcellnew++) 
        {
        cellarray->GetNextCell(npts,v_id);
        newCellId = tets->InsertNextCell(npts,v_id);
        outCD->CopyData(inCD,cellId,newCellId);
        }
      cellarray->Delete();
    }
  arraytetra->Delete();
}
