[occt.git] / src / Voxel / Voxel_Selector.cxx

// Created on: 2008-07-30
// Created by: Vladislav ROMASHKO
// Copyright (c) 2008-2014 OPEN CASCADE SAS
//
// This file is part of Open CASCADE Technology software library.
//
// This library is free software; you can redistribute it and/or modify it under
// the terms of the GNU Lesser General Public License version 2.1 as published
// by the Free Software Foundation, with special exception defined in the file
// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
// distribution for complete text of the license and disclaimer of any warranty.
//
// Alternatively, this file may be used under the terms of Open CASCADE
// commercial license or contractual agreement.

#include <Voxel_Selector.ixx>

#include <gp_Pnt.hxx>
#include <gp_Lin.hxx>
#include <gp_Dir.hxx>
#include <gp_Pln.hxx>

#include <IntAna_Quadric.hxx>
#include <IntAna_IntConicQuad.hxx>

#include <TColStd_MapOfInteger.hxx>

Voxel_Selector::Voxel_Selector():myVoxels(0)
{

}

Voxel_Selector::Voxel_Selector(const Handle(V3d_View)& view):myView(view),myVoxels(0)
{

}

void Voxel_Selector::Init(const Handle(V3d_View)& view)
{
  myView = view;
}

void Voxel_Selector::SetVoxels(const Voxel_BoolDS& voxels)
{
  myIsBool = 1;
  myVoxels = (void*) &voxels;
}

void Voxel_Selector::SetVoxels(const Voxel_ColorDS& voxels)
{
  myIsBool = 0;
  myVoxels = (void*) &voxels;
}

void Voxel_Selector::SetVoxels(const Voxel_ROctBoolDS& voxels)
{
  myIsBool = 2;
  myVoxels = (void*) &voxels;
}

// This function is copied from ViewerTest_RelationCommands.cxx
static Standard_Boolean ComputeIntersection(const gp_Lin& L,const gp_Pln& ThePl, gp_Pnt& TheInter)
{
  static IntAna_Quadric TheQuad;
  TheQuad.SetQuadric(ThePl);
  static IntAna_IntConicQuad QQ;
  QQ.Perform(L,TheQuad);
  if(QQ.IsDone()){
    if(QQ.NbPoints()>0){
      TheInter = QQ.Point(1);
      return Standard_True;
    }
  }
  return Standard_False;
}

static inline Standard_Integer GetIVoxel(const Standard_Integer ix, const Standard_Integer iy, const Standard_Integer iz, 
					 const Standard_Integer nbx, const Standard_Integer nbxy)
{
  return ix + iy * nbx + iz * nbxy;
}

static inline Standard_Boolean Get(const Standard_Address voxels, const Standard_Integer isBool,
				   const Standard_Integer ix, const Standard_Integer iy, const Standard_Integer iz)
{
  switch (isBool)
  {
    case 0:
      return ((Voxel_ColorDS*) voxels)->Get(ix, iy, iz) > 0;
    case 1:
      return ((Voxel_BoolDS*) voxels)->Get(ix, iy, iz);
    case 2:
    {
      Standard_Integer deepness = ((Voxel_ROctBoolDS*) voxels)->Deepness(ix, iy, iz);
      switch (deepness)
      {
	case 0:
	  return ((Voxel_ROctBoolDS*) voxels)->Get(ix, iy, iz);
	case 1:
	{
	  for (Standard_Integer i = 0; i < 8; i++)
	  {
	    if (((Voxel_ROctBoolDS*) voxels)->Get(ix, iy, iz, i) == Standard_True)
	      return Standard_True;
	  }
	  break;
	}
	case 2:
	{
	  for (Standard_Integer i = 0; i < 8; i++)
	  {
	    for (Standard_Integer j = 0; j < 8; j++)
	    {
	      if (((Voxel_ROctBoolDS*) voxels)->Get(ix, iy, iz, i, j) == Standard_True)
		return Standard_True;
	    }
	  }
	  break;
	}
      }
    }
  }
  return Standard_False;
}

Standard_Boolean Voxel_Selector::Detect(const Standard_Integer winx, const Standard_Integer winy,
					Standard_Integer& ixdetect, Standard_Integer& iydetect, Standard_Integer& izdetect)
{
  ixdetect = -1; iydetect = -1; izdetect = -1;
  if (myView.IsNull() || !myVoxels)
    return Standard_False;

  Voxel_DS* ds = 0;
  switch (myIsBool)
  {
    case 0:
      ds = (Voxel_ColorDS*) myVoxels;
      break;
    case 1:
      ds = (Voxel_BoolDS*) myVoxels;
      break;
    case 2:
      ds = (Voxel_ROctBoolDS*) myVoxels;
      break;
  }
  Standard_Integer nbx = ds->GetNbX(), nby = ds->GetNbY(), nbz = ds->GetNbZ(), nbxy = nbx * nby;

  // Construct a line perpendicular to the screen
  Standard_Real eyex, eyey, eyez, nx, ny, nz;
  myView->Convert(winx, winy, eyex, eyey, eyez);
  myView->Proj(nx, ny, nz);
  gp_Pnt peye(eyex, eyey, eyez);
  gp_Lin line(peye, gp_Dir(nx, ny, nz));

  // Find the first voxel meeting the line at entrance to the cube of voxels.
  // Construct planes of the cube of voxels
  Standard_Real xstart = ds->GetX(), ystart = ds->GetY(), zstart = ds->GetZ();
  Standard_Real xlen = ds->GetXLen(), ylen = ds->GetYLen(), zlen = ds->GetZLen();
  Standard_Real xend = xstart + xlen, yend = ystart + ylen, zend = zstart + zlen;
  gp_Pln xplane_minus(gp_Pnt(xstart, ystart, zstart), -gp::DX());
  gp_Pln xplane_plus (gp_Pnt(xend,   ystart, zstart),  gp::DX());
  gp_Pln yplane_minus(gp_Pnt(xstart, ystart, zstart), -gp::DY());
  gp_Pln yplane_plus (gp_Pnt(xstart, yend,   zstart),  gp::DY());
  gp_Pln zplane_minus(gp_Pnt(xstart, ystart, zstart), -gp::DZ());
  gp_Pln zplane_plus (gp_Pnt(xstart, ystart, zend),    gp::DZ());
  // Intersect the planes with the line.
  gp_Pnt pintersection, p;
  Standard_Real depth = DBL_MAX, d;
  Standard_Integer iplane = -1; // not found
  if (ComputeIntersection(line, xplane_minus, p)) // -X
  {
    if (p.Y() >= ystart && p.Y() <= yend &&
	p.Z() >= zstart && p.Z() <= zend)
    {
      p.SetX(xstart);
      depth = peye.SquareDistance(p);
      iplane = 0;
      pintersection = p;
    }
  }
  if (ComputeIntersection(line, xplane_plus, p)) // +X
  {
    if (p.Y() >= ystart && p.Y() <= yend &&
	p.Z() >= zstart && p.Z() <= zend)
    {
      d = peye.SquareDistance(p);
      if (d < depth)
      {
	p.SetX(xend);
	depth = d;
	iplane = 1;
	pintersection = p;
      }
    }
  }
  if (ComputeIntersection(line, yplane_minus, p)) // -Y
  {
    if (p.X() >= xstart && p.X() <= xend &&
	p.Z() >= zstart && p.Z() <= zend)
    {
      d = peye.SquareDistance(p);
      if (d < depth)
      {
	p.SetY(ystart);
	depth = d;
	iplane = 2;
	pintersection = p;
      }
    }
  }
  if (ComputeIntersection(line, yplane_plus, p)) // +Y
  {
    if (p.X() >= xstart && p.X() <= xend &&
	p.Z() >= zstart && p.Z() <= zend)
    {
      d = peye.SquareDistance(p);
      if (d < depth)
      {
	p.SetY(yend);
	depth = d;
	iplane = 3;
	pintersection = p;
      }
    }
  }
  if (ComputeIntersection(line, zplane_minus, p)) // -Z
  {
    if (p.X() >= xstart && p.X() <= xend &&
	p.Y() >= ystart && p.Y() <= yend)
    {
      d = peye.SquareDistance(p);
      if (d < depth)
      {
	p.SetZ(zstart);
	depth = d;
	iplane = 4;
	pintersection = p;
      }
    }
  }
  if (ComputeIntersection(line, zplane_plus, p)) // +Z
  {
    if (p.X() >= xstart && p.X() <= xend &&
	p.Y() >= ystart && p.Y() <= yend)
    {
      d = peye.SquareDistance(p);
      if (d < depth)
      {
	p.SetZ(zend);
	depth = d;
	iplane = 5;
	pintersection = p;
      }
    }
  }
  // Find the voxel on the detected plane
  if (iplane == -1)
    return Standard_False;
  Standard_Integer ix, iy, iz;
  if (!ds->GetVoxel(pintersection.X(), pintersection.Y(), pintersection.Z(), ix, iy, iz))
    return Standard_False;
  ixdetect = ix; iydetect = iy; izdetect = iz;

  // Find a non-zero voxel at the line
  Standard_Real xmin = xlen / Standard_Real(nbx),
                ymin = ylen / Standard_Real(nby),
                zmin = zlen / Standard_Real(nbz),
                vmin = sqrt(xmin * xmin + ymin * ymin + zmin * zmin) / 2.0;
  Standard_Real xc, yc, zc, dist, distmin = DBL_MAX;
  TColStd_MapOfInteger passed;
  while (!Get(myVoxels, myIsBool, ixdetect, iydetect, izdetect))
  {
    // Memorize already checked voxels
    if (!passed.Add(GetIVoxel(ixdetect, iydetect, izdetect, nbx, nbxy)))
      return Standard_False;

    distmin = DBL_MAX;
    ix = ixdetect; iy = iydetect; iz = izdetect;

    //1: -X neighbour
    if (ix - 1 >= 0 && !passed.Contains(GetIVoxel(ix - 1, iy, iz, nbx, nbxy)))
    {
      ds->GetCenter(ix - 1, iy, iz, xc, yc, zc);
      dist = line.Distance(gp_Pnt(xc, yc, zc));
      if (dist < vmin && dist < distmin)
      {
	ixdetect = ix - 1; iydetect = iy; izdetect = iz;
	distmin = dist;
      }
    }
    //2: +X neighbour
    if (ix + 1 < nbx && !passed.Contains(GetIVoxel(ix + 1, iy, iz, nbx, nbxy)))
    {
      ds->GetCenter(ix + 1, iy, iz, xc, yc, zc);
      dist = line.Distance(gp_Pnt(xc, yc, zc));
      if (dist < vmin && dist < distmin)
      {
	ixdetect = ix + 1; iydetect = iy; izdetect = iz;
	distmin = dist;
      }
    }
    //3: -Y neighbour
    if (iy - 1 >= 0 && !passed.Contains(GetIVoxel(ix, iy - 1, iz, nbx, nbxy)))
    {
      ds->GetCenter(ix, iy - 1, iz, xc, yc, zc);
      dist = line.Distance(gp_Pnt(xc, yc, zc));
      if (dist < vmin && dist < distmin)
      {
	ixdetect = ix; iydetect = iy - 1; izdetect = iz;
	distmin = dist;
      }
    }
    //4: +Y neighbour
    if (iy + 1 < nby && !passed.Contains(GetIVoxel(ix, iy + 1, iz, nbx, nbxy)))
    {
      ds->GetCenter(ix, iy + 1, iz, xc, yc, zc);
      dist = line.Distance(gp_Pnt(xc, yc, zc));
      if (dist < vmin && dist < distmin)
      {
	ixdetect = ix; iydetect = iy + 1; izdetect = iz;
	distmin = dist;
      }
    }
    //5: -Z neighbour
    if (iz - 1 >= 0 && !passed.Contains(GetIVoxel(ix, iy, iz - 1, nbx, nbxy)))
    {
      ds->GetCenter(ix, iy, iz - 1, xc, yc, zc);
      dist = line.Distance(gp_Pnt(xc, yc, zc));
      if (dist < vmin && dist < distmin)
      {
	ixdetect = ix; iydetect = iy; izdetect = iz - 1;
	distmin = dist;
      }
    }
    //6: +Z neighbour
    if (iz + 1 < nbz && !passed.Contains(GetIVoxel(ix, iy, iz + 1, nbx, nbxy)))
    {
      ds->GetCenter(ix, iy, iz + 1, xc, yc, zc);
      dist = line.Distance(gp_Pnt(xc, yc, zc));
      if (dist < vmin && dist < distmin)
      {
	ixdetect = ix; iydetect = iy; izdetect = iz + 1;
	distmin = dist;
      }
    }

    // Diagonal voxels
    //7: -X-Y neighbour
    if (ix - 1 >= 0 && iy - 1 >= 0 && !passed.Contains(GetIVoxel(ix - 1, iy - 1, iz, nbx, nbxy)))
    {
      ds->GetCenter(ix - 1, iy - 1, iz, xc, yc, zc);
      dist = line.Distance(gp_Pnt(xc, yc, zc));
      if (dist < vmin && dist < distmin)
      {
	ixdetect = ix - 1; iydetect = iy - 1; izdetect = iz;
	distmin = dist;
      }
    }
    //8: -X-Z neighbour
    if (ix - 1 >= 0 && iz - 1 >= 0 && !passed.Contains(GetIVoxel(ix - 1, iy, iz - 1, nbx, nbxy)))
    {
      ds->GetCenter(ix - 1, iy, iz - 1, xc, yc, zc);
      dist = line.Distance(gp_Pnt(xc, yc, zc));
      if (dist < vmin && dist < distmin)
      {
	ixdetect = ix - 1; iydetect = iy; izdetect = iz - 1;
	distmin = dist;
      }
    }
    //9: -Y-Z neighbour
    if (iy - 1 >= 0 && iz - 1 >= 0 && !passed.Contains(GetIVoxel(ix, iy - 1, iz - 1, nbx, nbxy)))
    {
      ds->GetCenter(ix, iy - 1, iz - 1, xc, yc, zc);
      dist = line.Distance(gp_Pnt(xc, yc, zc));
      if (dist < vmin && dist < distmin)
      {
	ixdetect = ix; iydetect = iy - 1; izdetect = iz - 1;
	distmin = dist;
      }
    }

    //10: +X-Y neighbour
    if (ix + 1 < nbx && iy - 1 >= 0 && !passed.Contains(GetIVoxel(ix + 1, iy - 1, iz, nbx, nbxy)))
    {
      ds->GetCenter(ix + 1, iy - 1, iz, xc, yc, zc);
      dist = line.Distance(gp_Pnt(xc, yc, zc));
      if (dist < vmin && dist < distmin)
      {
	ixdetect = ix + 1; iydetect = iy - 1; izdetect = iz;
	distmin = dist;
      }
    }
    //11: +X-Z neighbour
    if (ix + 1 < nbx && iz - 1 >= 0 && !passed.Contains(GetIVoxel(ix + 1, iy, iz - 1, nbx, nbxy)))
    {
      ds->GetCenter(ix + 1, iy, iz - 1, xc, yc, zc);
      dist = line.Distance(gp_Pnt(xc, yc, zc));
      if (dist < vmin && dist < distmin)
      {
	ixdetect = ix + 1; iydetect = iy; izdetect = iz - 1;
	distmin = dist;
      }
    }
    //12: +Y-Z neighbour
    if (iy + 1 < nby && iz - 1 >= 0 && !passed.Contains(GetIVoxel(ix, iy + 1, iz - 1, nbx, nbxy)))
    {
      ds->GetCenter(ix, iy + 1, iz - 1, xc, yc, zc);
      dist = line.Distance(gp_Pnt(xc, yc, zc));
      if (dist < vmin && dist < distmin)
      {
	ixdetect = ix; iydetect = iy + 1; izdetect = iz - 1;
	distmin = dist;
      }
    }

    //13: -X+Y neighbour
    if (ix - 1 >= 0 && iy + 1 < nby && !passed.Contains(GetIVoxel(ix - 1, iy + 1, iz, nbx, nbxy)))
    {
      ds->GetCenter(ix - 1, iy + 1, iz, xc, yc, zc);
      dist = line.Distance(gp_Pnt(xc, yc, zc));
      if (dist < vmin && dist < distmin)
      {
	ixdetect = ix - 1; iydetect = iy + 1; izdetect = iz;
	distmin = dist;
      }
    }
    //14: -X+Z neighbour
    if (ix - 1 >= 0 && iz + 1 < nbz && !passed.Contains(GetIVoxel(ix - 1, iy, iz + 1, nbx, nbxy)))
    {
      ds->GetCenter(ix - 1, iy, iz + 1, xc, yc, zc);
      dist = line.Distance(gp_Pnt(xc, yc, zc));
      if (dist < vmin && dist < distmin)
      {
	ixdetect = ix - 1; iydetect = iy; izdetect = iz + 1;
	distmin = dist;
      }
    }
    //15: -Y+Z neighbour
    if (iy - 1 >= 0 && iz + 1 < nbz && !passed.Contains(GetIVoxel(ix, iy - 1, iz + 1, nbx, nbxy)))
    {
      ds->GetCenter(ix, iy - 1, iz + 1, xc, yc, zc);
      dist = line.Distance(gp_Pnt(xc, yc, zc));
      if (dist < vmin && dist < distmin)
      {
	ixdetect = ix; iydetect = iy - 1; izdetect = iz + 1;
	distmin = dist;
      }
    }

    //16: +X+Y neighbour
    if (ix + 1 < nbx && iy + 1 < nby && !passed.Contains(GetIVoxel(ix + 1, iy + 1, iz, nbx, nbxy)))
    {
      ds->GetCenter(ix + 1, iy + 1, iz, xc, yc, zc);
      dist = line.Distance(gp_Pnt(xc, yc, zc));
      if (dist < vmin && dist < distmin)
      {
	ixdetect = ix + 1; iydetect = iy + 1; izdetect = iz;
	distmin = dist;
      }
    }
    //17: +X+Z neighbour
    if (ix + 1 < nbx && iz + 1 < nbz && !passed.Contains(GetIVoxel(ix + 1, iy, iz + 1, nbx, nbxy)))
    {
      ds->GetCenter(ix + 1, iy, iz + 1, xc, yc, zc);
      dist = line.Distance(gp_Pnt(xc, yc, zc));
      if (dist < vmin && dist < distmin)
      {
	ixdetect = ix + 1; iydetect = iy; izdetect = iz + 1;
	distmin = dist;
      }
    }
    //18: +Y+Z neighbour
    if (iy + 1 < nby && iz + 1 < nbz && !passed.Contains(GetIVoxel(ix, iy + 1, iz + 1, nbx, nbxy)))
    {
      ds->GetCenter(ix, iy + 1, iz + 1, xc, yc, zc);
      dist = line.Distance(gp_Pnt(xc, yc, zc));
      if (dist < vmin && dist < distmin)
      {
	ixdetect = ix; iydetect = iy + 1; izdetect = iz + 1;
	distmin = dist;
      }
    }

    // Farest neighbours
    //19: -X-Y-Z neighbour
    if (ix - 1 >= 0 && iy - 1 >= 0 && iz - 1 >= 0 && !passed.Contains(GetIVoxel(ix - 1, iy - 1, iz - 1, nbx, nbxy)))
    {
      ds->GetCenter(ix - 1, iy - 1, iz - 1, xc, yc, zc);
      dist = line.Distance(gp_Pnt(xc, yc, zc));
      if (dist < vmin && dist < distmin)
      {
	ixdetect = ix - 1; iydetect = iy - 1; izdetect = iz - 1;
	distmin = dist;
      }
    }
    //20: +X-Y-Z neighbour
    if (ix + 1 < nbx && iy - 1 >= 0 && iz - 1 >= 0 && !passed.Contains(GetIVoxel(ix + 1, iy - 1, iz - 1, nbx, nbxy)))
    {
      ds->GetCenter(ix + 1, iy - 1, iz - 1, xc, yc, zc);
      dist = line.Distance(gp_Pnt(xc, yc, zc));
      if (dist < vmin && dist < distmin)
      {
	ixdetect = ix + 1; iydetect = iy - 1; izdetect = iz - 1;
 	distmin = dist;
     }
    }
    //21: -X+Y-Z neighbour
    if (ix - 1 >= 0 && iy + 1 < nby && iz - 1 >= 0 && !passed.Contains(GetIVoxel(ix - 1, iy + 1, iz - 1, nbx, nbxy)))
    {
      ds->GetCenter(ix - 1, iy + 1, iz - 1, xc, yc, zc);
      dist = line.Distance(gp_Pnt(xc, yc, zc));
      if (dist < vmin && dist < distmin)
      {
	ixdetect = ix - 1; iydetect = iy + 1; izdetect = iz - 1;
	distmin = dist;
      }
    }
    //22: -X-Y+Z neighbour
    if (ix - 1 >= 0 && iy - 1 >= 0 && iz + 1 < nbz && !passed.Contains(GetIVoxel(ix - 1, iy - 1, iz + 1, nbx, nbxy)))
    {
      ds->GetCenter(ix - 1, iy - 1, iz + 1, xc, yc, zc);
      dist = line.Distance(gp_Pnt(xc, yc, zc));
      if (dist < vmin && dist < distmin)
      {
	ixdetect = ix - 1; iydetect = iy - 1; izdetect = iz + 1;
	distmin = dist;
      }
    }
    //23: +X+Y-Z neighbour
    if (ix + 1 < nbx && iy + 1 < nby && iz - 1 >= 0 && !passed.Contains(GetIVoxel(ix + 1, iy + 1, iz - 1, nbx, nbxy)))
    {
      ds->GetCenter(ix + 1, iy + 1, iz - 1, xc, yc, zc);
      dist = line.Distance(gp_Pnt(xc, yc, zc));
      if (dist < vmin && dist < distmin)
      {
	ixdetect = ix + 1; iydetect = iy + 1; izdetect = iz - 1;
	distmin = dist;
      }
    }
    //24: +X-Y+Z neighbour
    if (ix + 1 < nbx && iy - 1 >= 0 && iz + 1 < nbz && !passed.Contains(GetIVoxel(ix + 1, iy - 1, iz + 1, nbx, nbxy)))
    {
      ds->GetCenter(ix + 1, iy - 1, iz + 1, xc, yc, zc);
      dist = line.Distance(gp_Pnt(xc, yc, zc));
      if (dist < vmin && dist < distmin)
      {
	ixdetect = ix + 1; iydetect = iy - 1; izdetect = iz + 1;
	distmin = dist;
      }
    }
    //25: -X+Y+Z neighbour
    if (ix - 1 >= 0 && iy + 1 < nby && iz + 1 < nbz && !passed.Contains(GetIVoxel(ix - 1, iy + 1, iz + 1, nbx, nbxy)))
    {
      ds->GetCenter(ix - 1, iy + 1, iz + 1, xc, yc, zc);
      dist = line.Distance(gp_Pnt(xc, yc, zc));
      if (dist < vmin && dist < distmin)
      {
	ixdetect = ix - 1; iydetect = iy + 1; izdetect = iz + 1;
	distmin = dist;
      }
    }
    //26: +X+Y+Z neighbour
    if (ix + 1 < nbx && iy + 1 < nby && iz + 1 < nbz && !passed.Contains(GetIVoxel(ix + 1, iy + 1, iz + 1, nbx, nbxy)))
    {
      ds->GetCenter(ix + 1, iy + 1, iz + 1, xc, yc, zc);
      dist = line.Distance(gp_Pnt(xc, yc, zc));
      if (dist < vmin && dist < distmin)
      {
	ixdetect = ix + 1; iydetect = iy + 1; izdetect = iz + 1;
	distmin = dist;
      }
    }

  } // End of while (zero-voxel...

  if (!Get(myVoxels, myIsBool, ixdetect, iydetect, izdetect))
    return Standard_False;
  return Standard_True;
}
Commit	Line	Data
b311480e	1	// Created on: 2008-07-30
b311480e	2	// Created by: Vladislav ROMASHKO
973c2be1	3	// Copyright (c) 2008-2014 OPEN CASCADE SAS
b311480e	4	//
973c2be1	5	// This file is part of Open CASCADE Technology software library.
b311480e	6	//
d5f74e42	7	// This library is free software; you can redistribute it and/or modify it under
d5f74e42	8	// the terms of the GNU Lesser General Public License version 2.1 as published
973c2be1	9	// by the Free Software Foundation, with special exception defined in the file
	10	// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
	11	// distribution for complete text of the license and disclaimer of any warranty.
b311480e	12	//
973c2be1	13	// Alternatively, this file may be used under the terms of Open CASCADE
973c2be1	14	// commercial license or contractual agreement.
7fd59977	15
	16	#include <Voxel_Selector.ixx>
	17
	18	#include <gp_Pnt.hxx>
	19	#include <gp_Lin.hxx>
	20	#include <gp_Dir.hxx>
	21	#include <gp_Pln.hxx>
	22
	23	#include <IntAna_Quadric.hxx>
	24	#include <IntAna_IntConicQuad.hxx>
	25
	26	#include <TColStd_MapOfInteger.hxx>
	27
	28	Voxel_Selector::Voxel_Selector():myVoxels(0)
	29	{
	30
	31	}
	32
	33	Voxel_Selector::Voxel_Selector(const Handle(V3d_View)& view):myView(view),myVoxels(0)
	34	{
	35
	36	}
	37
	38	void Voxel_Selector::Init(const Handle(V3d_View)& view)
	39	{
	40	myView = view;
	41	}
	42
	43	void Voxel_Selector::SetVoxels(const Voxel_BoolDS& voxels)
	44	{
	45	myIsBool = 1;
	46	myVoxels = (void*) &voxels;
	47	}
	48
	49	void Voxel_Selector::SetVoxels(const Voxel_ColorDS& voxels)
	50	{
	51	myIsBool = 0;
	52	myVoxels = (void*) &voxels;
	53	}
	54
	55	void Voxel_Selector::SetVoxels(const Voxel_ROctBoolDS& voxels)
	56	{
	57	myIsBool = 2;
	58	myVoxels = (void*) &voxels;
	59	}
	60
	61	// This function is copied from ViewerTest_RelationCommands.cxx
	62	static Standard_Boolean ComputeIntersection(const gp_Lin& L,const gp_Pln& ThePl, gp_Pnt& TheInter)
	63	{
	64	static IntAna_Quadric TheQuad;
	65	TheQuad.SetQuadric(ThePl);
	66	static IntAna_IntConicQuad QQ;
	67	QQ.Perform(L,TheQuad);
	68	if(QQ.IsDone()){
	69	if(QQ.NbPoints()>0){
	70	TheInter = QQ.Point(1);
	71	return Standard_True;
	72	}
	73	}
	74	return Standard_False;
	75	}
	76
	77	static inline Standard_Integer GetIVoxel(const Standard_Integer ix, const Standard_Integer iy, const Standard_Integer iz,
	78	const Standard_Integer nbx, const Standard_Integer nbxy)
79	{
80	return ix + iy * nbx + iz * nbxy;
81	}
82
83	static inline Standard_Boolean Get(const Standard_Address voxels, const Standard_Integer isBool,
84	const Standard_Integer ix, const Standard_Integer iy, const Standard_Integer iz)
85	{
86	switch (isBool)
87	{
88	case 0:
89	return ((Voxel_ColorDS*) voxels)->Get(ix, iy, iz) > 0;
90	case 1:
91	return ((Voxel_BoolDS*) voxels)->Get(ix, iy, iz);
92	case 2:
93	{
94	Standard_Integer deepness = ((Voxel_ROctBoolDS*) voxels)->Deepness(ix, iy, iz);
95	switch (deepness)
96	{
97	case 0:
98	return ((Voxel_ROctBoolDS*) voxels)->Get(ix, iy, iz);
99	case 1:
100	{
101	for (Standard_Integer i = 0; i < 8; i++)
102	{
103	if (((Voxel_ROctBoolDS*) voxels)->Get(ix, iy, iz, i) == Standard_True)
104	return Standard_True;
105	}
106	break;
107	}
108	case 2:
109	{
110	for (Standard_Integer i = 0; i < 8; i++)
111	{
112	for (Standard_Integer j = 0; j < 8; j++)
113	{
114	if (((Voxel_ROctBoolDS*) voxels)->Get(ix, iy, iz, i, j) == Standard_True)
115	return Standard_True;
116	}
117	}
118	break;
119	}
120	}
121	}
122	}
123	return Standard_False;
124	}
125
126	Standard_Boolean Voxel_Selector::Detect(const Standard_Integer winx, const Standard_Integer winy,
127	Standard_Integer& ixdetect, Standard_Integer& iydetect, Standard_Integer& izdetect)
128	{
129	ixdetect = -1; iydetect = -1; izdetect = -1;
130	if (myView.IsNull() \|\| !myVoxels)
131	return Standard_False;
132
133	Voxel_DS* ds = 0;
134	switch (myIsBool)
135	{
136	case 0:
137	ds = (Voxel_ColorDS*) myVoxels;
138	break;
139	case 1:
140	ds = (Voxel_BoolDS*) myVoxels;
141	break;
142	case 2:
143	ds = (Voxel_ROctBoolDS*) myVoxels;
144	break;
145	}
146	Standard_Integer nbx = ds->GetNbX(), nby = ds->GetNbY(), nbz = ds->GetNbZ(), nbxy = nbx * nby;
147
148	// Construct a line perpendicular to the screen
149	Standard_Real eyex, eyey, eyez, nx, ny, nz;
150	myView->Convert(winx, winy, eyex, eyey, eyez);
151	myView->Proj(nx, ny, nz);
152	gp_Pnt peye(eyex, eyey, eyez);
153	gp_Lin line(peye, gp_Dir(nx, ny, nz));
154
155	// Find the first voxel meeting the line at entrance to the cube of voxels.
156	// Construct planes of the cube of voxels
157	Standard_Real xstart = ds->GetX(), ystart = ds->GetY(), zstart = ds->GetZ();
158	Standard_Real xlen = ds->GetXLen(), ylen = ds->GetYLen(), zlen = ds->GetZLen();
159	Standard_Real xend = xstart + xlen, yend = ystart + ylen, zend = zstart + zlen;
160	gp_Pln xplane_minus(gp_Pnt(xstart, ystart, zstart), -gp::DX());
161	gp_Pln xplane_plus (gp_Pnt(xend, ystart, zstart), gp::DX());
162	gp_Pln yplane_minus(gp_Pnt(xstart, ystart, zstart), -gp::DY());
163	gp_Pln yplane_plus (gp_Pnt(xstart, yend, zstart), gp::DY());
164	gp_Pln zplane_minus(gp_Pnt(xstart, ystart, zstart), -gp::DZ());
165	gp_Pln zplane_plus (gp_Pnt(xstart, ystart, zend), gp::DZ());
166	// Intersect the planes with the line.
167	gp_Pnt pintersection, p;
168	Standard_Real depth = DBL_MAX, d;
169	Standard_Integer iplane = -1; // not found
170	if (ComputeIntersection(line, xplane_minus, p)) // -X
171	{
172	if (p.Y() >= ystart && p.Y() <= yend &&
173	p.Z() >= zstart && p.Z() <= zend)
174	{
175	p.SetX(xstart);
176	depth = peye.SquareDistance(p);
177	iplane = 0;
178	pintersection = p;
179	}
180	}
181	if (ComputeIntersection(line, xplane_plus, p)) // +X
182	{
183	if (p.Y() >= ystart && p.Y() <= yend &&
184	p.Z() >= zstart && p.Z() <= zend)
185	{
186	d = peye.SquareDistance(p);
187	if (d < depth)
188	{
189	p.SetX(xend);
190	depth = d;
191	iplane = 1;
192	pintersection = p;
193	}
194	}
195	}
196	if (ComputeIntersection(line, yplane_minus, p)) // -Y
197	{
198	if (p.X() >= xstart && p.X() <= xend &&
199	p.Z() >= zstart && p.Z() <= zend)
200	{
201	d = peye.SquareDistance(p);
202	if (d < depth)
203	{
204	p.SetY(ystart);
205	depth = d;
206	iplane = 2;
207	pintersection = p;
208	}
209	}
210	}
211	if (ComputeIntersection(line, yplane_plus, p)) // +Y
212	{
213	if (p.X() >= xstart && p.X() <= xend &&
214	p.Z() >= zstart && p.Z() <= zend)
215	{
216	d = peye.SquareDistance(p);
217	if (d < depth)
218	{
219	p.SetY(yend);
220	depth = d;
221	iplane = 3;
222	pintersection = p;
223	}
224	}
225	}
226	if (ComputeIntersection(line, zplane_minus, p)) // -Z
227	{
228	if (p.X() >= xstart && p.X() <= xend &&
229	p.Y() >= ystart && p.Y() <= yend)
230	{
231	d = peye.SquareDistance(p);
232	if (d < depth)
233	{
234	p.SetZ(zstart);
235	depth = d;
236	iplane = 4;
237	pintersection = p;
238	}
239	}
240	}
241	if (ComputeIntersection(line, zplane_plus, p)) // +Z
242	{
243	if (p.X() >= xstart && p.X() <= xend &&
244	p.Y() >= ystart && p.Y() <= yend)
245	{
246	d = peye.SquareDistance(p);
247	if (d < depth)
248	{
249	p.SetZ(zend);
250	depth = d;
251	iplane = 5;
252	pintersection = p;
253	}
254	}
255	}
256	// Find the voxel on the detected plane
257	if (iplane == -1)
258	return Standard_False;
259	Standard_Integer ix, iy, iz;
260	if (!ds->GetVoxel(pintersection.X(), pintersection.Y(), pintersection.Z(), ix, iy, iz))
261	return Standard_False;
262	ixdetect = ix; iydetect = iy; izdetect = iz;
263
264	// Find a non-zero voxel at the line
265	Standard_Real xmin = xlen / Standard_Real(nbx),
266	ymin = ylen / Standard_Real(nby),
267	zmin = zlen / Standard_Real(nbz),
268	vmin = sqrt(xmin * xmin + ymin * ymin + zmin * zmin) / 2.0;
269	Standard_Real xc, yc, zc, dist, distmin = DBL_MAX;
270	TColStd_MapOfInteger passed;
271	while (!Get(myVoxels, myIsBool, ixdetect, iydetect, izdetect))
272	{
273	// Memorize already checked voxels
274	if (!passed.Add(GetIVoxel(ixdetect, iydetect, izdetect, nbx, nbxy)))
275	return Standard_False;
276
277	distmin = DBL_MAX;
278	ix = ixdetect; iy = iydetect; iz = izdetect;
279
280	//1: -X neighbour
281	if (ix - 1 >= 0 && !passed.Contains(GetIVoxel(ix - 1, iy, iz, nbx, nbxy)))
282	{
283	ds->GetCenter(ix - 1, iy, iz, xc, yc, zc);
284	dist = line.Distance(gp_Pnt(xc, yc, zc));
285	if (dist < vmin && dist < distmin)
286	{
287	ixdetect = ix - 1; iydetect = iy; izdetect = iz;
288	distmin = dist;
289	}
290	}
291	//2: +X neighbour
292	if (ix + 1 < nbx && !passed.Contains(GetIVoxel(ix + 1, iy, iz, nbx, nbxy)))
293	{
294	ds->GetCenter(ix + 1, iy, iz, xc, yc, zc);
295	dist = line.Distance(gp_Pnt(xc, yc, zc));
296	if (dist < vmin && dist < distmin)
297	{
298	ixdetect = ix + 1; iydetect = iy; izdetect = iz;
299	distmin = dist;
300	}
301	}
302	//3: -Y neighbour
303	if (iy - 1 >= 0 && !passed.Contains(GetIVoxel(ix, iy - 1, iz, nbx, nbxy)))
304	{
305	ds->GetCenter(ix, iy - 1, iz, xc, yc, zc);
306	dist = line.Distance(gp_Pnt(xc, yc, zc));
307	if (dist < vmin && dist < distmin)
308	{
309	ixdetect = ix; iydetect = iy - 1; izdetect = iz;
310	distmin = dist;
311	}
312	}
313	//4: +Y neighbour
314	if (iy + 1 < nby && !passed.Contains(GetIVoxel(ix, iy + 1, iz, nbx, nbxy)))
315	{
316	ds->GetCenter(ix, iy + 1, iz, xc, yc, zc);
317	dist = line.Distance(gp_Pnt(xc, yc, zc));
318	if (dist < vmin && dist < distmin)
319	{
320	ixdetect = ix; iydetect = iy + 1; izdetect = iz;
321	distmin = dist;
322	}
323	}
324	//5: -Z neighbour
325	if (iz - 1 >= 0 && !passed.Contains(GetIVoxel(ix, iy, iz - 1, nbx, nbxy)))
326	{
327	ds->GetCenter(ix, iy, iz - 1, xc, yc, zc);
328	dist = line.Distance(gp_Pnt(xc, yc, zc));
329	if (dist < vmin && dist < distmin)
330	{
331	ixdetect = ix; iydetect = iy; izdetect = iz - 1;
332	distmin = dist;
333	}
334	}
335	//6: +Z neighbour
336	if (iz + 1 < nbz && !passed.Contains(GetIVoxel(ix, iy, iz + 1, nbx, nbxy)))
337	{
338	ds->GetCenter(ix, iy, iz + 1, xc, yc, zc);
339	dist = line.Distance(gp_Pnt(xc, yc, zc));
340	if (dist < vmin && dist < distmin)
341	{
342	ixdetect = ix; iydetect = iy; izdetect = iz + 1;
343	distmin = dist;
344	}
345	}
346
347	// Diagonal voxels
348	//7: -X-Y neighbour
349	if (ix - 1 >= 0 && iy - 1 >= 0 && !passed.Contains(GetIVoxel(ix - 1, iy - 1, iz, nbx, nbxy)))
350	{
351	ds->GetCenter(ix - 1, iy - 1, iz, xc, yc, zc);
352	dist = line.Distance(gp_Pnt(xc, yc, zc));
353	if (dist < vmin && dist < distmin)
354	{
355	ixdetect = ix - 1; iydetect = iy - 1; izdetect = iz;
356	distmin = dist;
357	}
358	}
359	//8: -X-Z neighbour
360	if (ix - 1 >= 0 && iz - 1 >= 0 && !passed.Contains(GetIVoxel(ix - 1, iy, iz - 1, nbx, nbxy)))
361	{
362	ds->GetCenter(ix - 1, iy, iz - 1, xc, yc, zc);
363	dist = line.Distance(gp_Pnt(xc, yc, zc));
364	if (dist < vmin && dist < distmin)
365	{
366	ixdetect = ix - 1; iydetect = iy; izdetect = iz - 1;
367	distmin = dist;
368	}
369	}
370	//9: -Y-Z neighbour
371	if (iy - 1 >= 0 && iz - 1 >= 0 && !passed.Contains(GetIVoxel(ix, iy - 1, iz - 1, nbx, nbxy)))
372	{
373	ds->GetCenter(ix, iy - 1, iz - 1, xc, yc, zc);
374	dist = line.Distance(gp_Pnt(xc, yc, zc));
375	if (dist < vmin && dist < distmin)
376	{
377	ixdetect = ix; iydetect = iy - 1; izdetect = iz - 1;
378	distmin = dist;
379	}
380	}
381
382	//10: +X-Y neighbour
383	if (ix + 1 < nbx && iy - 1 >= 0 && !passed.Contains(GetIVoxel(ix + 1, iy - 1, iz, nbx, nbxy)))
384	{
385	ds->GetCenter(ix + 1, iy - 1, iz, xc, yc, zc);
386	dist = line.Distance(gp_Pnt(xc, yc, zc));
387	if (dist < vmin && dist < distmin)
388	{
389	ixdetect = ix + 1; iydetect = iy - 1; izdetect = iz;
390	distmin = dist;
391	}
392	}
393	//11: +X-Z neighbour
394	if (ix + 1 < nbx && iz - 1 >= 0 && !passed.Contains(GetIVoxel(ix + 1, iy, iz - 1, nbx, nbxy)))
395	{
396	ds->GetCenter(ix + 1, iy, iz - 1, xc, yc, zc);
397	dist = line.Distance(gp_Pnt(xc, yc, zc));
398	if (dist < vmin && dist < distmin)
399	{
400	ixdetect = ix + 1; iydetect = iy; izdetect = iz - 1;
401	distmin = dist;
402	}
403	}
404	//12: +Y-Z neighbour
405	if (iy + 1 < nby && iz - 1 >= 0 && !passed.Contains(GetIVoxel(ix, iy + 1, iz - 1, nbx, nbxy)))
406	{
407	ds->GetCenter(ix, iy + 1, iz - 1, xc, yc, zc);
408	dist = line.Distance(gp_Pnt(xc, yc, zc));
409	if (dist < vmin && dist < distmin)
410	{
411	ixdetect = ix; iydetect = iy + 1; izdetect = iz - 1;
412	distmin = dist;
413	}
414	}
415
416	//13: -X+Y neighbour
417	if (ix - 1 >= 0 && iy + 1 < nby && !passed.Contains(GetIVoxel(ix - 1, iy + 1, iz, nbx, nbxy)))
418	{
419	ds->GetCenter(ix - 1, iy + 1, iz, xc, yc, zc);
420	dist = line.Distance(gp_Pnt(xc, yc, zc));
421	if (dist < vmin && dist < distmin)
422	{
423	ixdetect = ix - 1; iydetect = iy + 1; izdetect = iz;
424	distmin = dist;
425	}
426	}
427	//14: -X+Z neighbour
428	if (ix - 1 >= 0 && iz + 1 < nbz && !passed.Contains(GetIVoxel(ix - 1, iy, iz + 1, nbx, nbxy)))
429	{
430	ds->GetCenter(ix - 1, iy, iz + 1, xc, yc, zc);
431	dist = line.Distance(gp_Pnt(xc, yc, zc));
432	if (dist < vmin && dist < distmin)
433	{
434	ixdetect = ix - 1; iydetect = iy; izdetect = iz + 1;
435	distmin = dist;
436	}
437	}
438	//15: -Y+Z neighbour
439	if (iy - 1 >= 0 && iz + 1 < nbz && !passed.Contains(GetIVoxel(ix, iy - 1, iz + 1, nbx, nbxy)))
440	{
441	ds->GetCenter(ix, iy - 1, iz + 1, xc, yc, zc);
442	dist = line.Distance(gp_Pnt(xc, yc, zc));
443	if (dist < vmin && dist < distmin)
444	{
445	ixdetect = ix; iydetect = iy - 1; izdetect = iz + 1;
446	distmin = dist;
447	}
448	}
449
450	//16: +X+Y neighbour
451	if (ix + 1 < nbx && iy + 1 < nby && !passed.Contains(GetIVoxel(ix + 1, iy + 1, iz, nbx, nbxy)))
452	{
453	ds->GetCenter(ix + 1, iy + 1, iz, xc, yc, zc);
454	dist = line.Distance(gp_Pnt(xc, yc, zc));
455	if (dist < vmin && dist < distmin)
456	{
457	ixdetect = ix + 1; iydetect = iy + 1; izdetect = iz;
458	distmin = dist;
459	}
460	}
461	//17: +X+Z neighbour
462	if (ix + 1 < nbx && iz + 1 < nbz && !passed.Contains(GetIVoxel(ix + 1, iy, iz + 1, nbx, nbxy)))
463	{
464	ds->GetCenter(ix + 1, iy, iz + 1, xc, yc, zc);
465	dist = line.Distance(gp_Pnt(xc, yc, zc));
466	if (dist < vmin && dist < distmin)
467	{
468	ixdetect = ix + 1; iydetect = iy; izdetect = iz + 1;
469	distmin = dist;
470	}
471	}
472	//18: +Y+Z neighbour
473	if (iy + 1 < nby && iz + 1 < nbz && !passed.Contains(GetIVoxel(ix, iy + 1, iz + 1, nbx, nbxy)))
474	{
475	ds->GetCenter(ix, iy + 1, iz + 1, xc, yc, zc);
476	dist = line.Distance(gp_Pnt(xc, yc, zc));
477	if (dist < vmin && dist < distmin)
478	{
479	ixdetect = ix; iydetect = iy + 1; izdetect = iz + 1;
480	distmin = dist;
481	}
482	}
483
484	// Farest neighbours
485	//19: -X-Y-Z neighbour
486	if (ix - 1 >= 0 && iy - 1 >= 0 && iz - 1 >= 0 && !passed.Contains(GetIVoxel(ix - 1, iy - 1, iz - 1, nbx, nbxy)))
487	{
488	ds->GetCenter(ix - 1, iy - 1, iz - 1, xc, yc, zc);
489	dist = line.Distance(gp_Pnt(xc, yc, zc));
490	if (dist < vmin && dist < distmin)
491	{
492	ixdetect = ix - 1; iydetect = iy - 1; izdetect = iz - 1;
493	distmin = dist;
494	}
495	}
496	//20: +X-Y-Z neighbour
497	if (ix + 1 < nbx && iy - 1 >= 0 && iz - 1 >= 0 && !passed.Contains(GetIVoxel(ix + 1, iy - 1, iz - 1, nbx, nbxy)))
498	{
499	ds->GetCenter(ix + 1, iy - 1, iz - 1, xc, yc, zc);
500	dist = line.Distance(gp_Pnt(xc, yc, zc));
501	if (dist < vmin && dist < distmin)
502	{
503	ixdetect = ix + 1; iydetect = iy - 1; izdetect = iz - 1;
504	distmin = dist;
505	}
506	}
507	//21: -X+Y-Z neighbour
508	if (ix - 1 >= 0 && iy + 1 < nby && iz - 1 >= 0 && !passed.Contains(GetIVoxel(ix - 1, iy + 1, iz - 1, nbx, nbxy)))
509	{
510	ds->GetCenter(ix - 1, iy + 1, iz - 1, xc, yc, zc);
511	dist = line.Distance(gp_Pnt(xc, yc, zc));
512	if (dist < vmin && dist < distmin)
513	{
514	ixdetect = ix - 1; iydetect = iy + 1; izdetect = iz - 1;
515	distmin = dist;
516	}
517	}
518	//22: -X-Y+Z neighbour
519	if (ix - 1 >= 0 && iy - 1 >= 0 && iz + 1 < nbz && !passed.Contains(GetIVoxel(ix - 1, iy - 1, iz + 1, nbx, nbxy)))
520	{
521	ds->GetCenter(ix - 1, iy - 1, iz + 1, xc, yc, zc);
522	dist = line.Distance(gp_Pnt(xc, yc, zc));
523	if (dist < vmin && dist < distmin)
524	{
525	ixdetect = ix - 1; iydetect = iy - 1; izdetect = iz + 1;
526	distmin = dist;
527	}
528	}
529	//23: +X+Y-Z neighbour
530	if (ix + 1 < nbx && iy + 1 < nby && iz - 1 >= 0 && !passed.Contains(GetIVoxel(ix + 1, iy + 1, iz - 1, nbx, nbxy)))
531	{
532	ds->GetCenter(ix + 1, iy + 1, iz - 1, xc, yc, zc);
533	dist = line.Distance(gp_Pnt(xc, yc, zc));
534	if (dist < vmin && dist < distmin)
535	{
536	ixdetect = ix + 1; iydetect = iy + 1; izdetect = iz - 1;
537	distmin = dist;
538	}
539	}
540	//24: +X-Y+Z neighbour
541	if (ix + 1 < nbx && iy - 1 >= 0 && iz + 1 < nbz && !passed.Contains(GetIVoxel(ix + 1, iy - 1, iz + 1, nbx, nbxy)))
542	{
543	ds->GetCenter(ix + 1, iy - 1, iz + 1, xc, yc, zc);
544	dist = line.Distance(gp_Pnt(xc, yc, zc));
545	if (dist < vmin && dist < distmin)
546	{
547	ixdetect = ix + 1; iydetect = iy - 1; izdetect = iz + 1;
548	distmin = dist;
549	}
550	}
551	//25: -X+Y+Z neighbour
552	if (ix - 1 >= 0 && iy + 1 < nby && iz + 1 < nbz && !passed.Contains(GetIVoxel(ix - 1, iy + 1, iz + 1, nbx, nbxy)))
553	{
554	ds->GetCenter(ix - 1, iy + 1, iz + 1, xc, yc, zc);
555	dist = line.Distance(gp_Pnt(xc, yc, zc));
556	if (dist < vmin && dist < distmin)
557	{
558	ixdetect = ix - 1; iydetect = iy + 1; izdetect = iz + 1;
559	distmin = dist;
560	}
561	}
562	//26: +X+Y+Z neighbour
563	if (ix + 1 < nbx && iy + 1 < nby && iz + 1 < nbz && !passed.Contains(GetIVoxel(ix + 1, iy + 1, iz + 1, nbx, nbxy)))
564	{
565	ds->GetCenter(ix + 1, iy + 1, iz + 1, xc, yc, zc);
566	dist = line.Distance(gp_Pnt(xc, yc, zc));
567	if (dist < vmin && dist < distmin)
568	{
569	ixdetect = ix + 1; iydetect = iy + 1; izdetect = iz + 1;
570	distmin = dist;
571	}
572	}
573
574	} // End of while (zero-voxel...
575
576	if (!Get(myVoxels, myIsBool, ixdetect, iydetect, izdetect))
577	return Standard_False;
578	return Standard_True;
579	}