[u/mrichter/AliRoot.git] / EVE / Reve / NLTProjector.cxx

#include "NLTProjector.h"
#include "RGTopFrame.h"
#include "NLTPolygonSet.h"
#include "PODs.h"
#include "PointSet.h"
#include "Track.h"

#include "TBuffer3D.h"
#include "TColor.h"
#include "TPointSet3D.h"
#include "TGeoManager.h"
#include "TGeoMatrix.h"

#include <list>

using namespace Reve;

namespace {
  struct Seg {
    Int_t v1;
    Int_t v2;

    Seg(Int_t i1=-1, Int_t i2=-1):v1(i1), v2(i2){};
  };
   
  typedef std::list<Seg>::iterator It_t;    
}

//______________________________________________________________________________
Vector*  NLTProjection::Project(Vector* origPnts, Int_t Npnts, Bool_t copy)
{
  Vector* pnts = 0; 
  if(copy) 
  {
    pnts = (Vector* )malloc(Npnts*sizeof(Vector));
    memcpy(pnts, origPnts, Npnts*sizeof(Vector));
  }
  else
  { 
    pnts =  origPnts;
  }
  return pnts;
}

//______________________________________________________________________________
Vector*  RhoZ::Project(Vector* origPnts, Int_t Npnts, Bool_t copy)
{
  Vector* pnts = NLTProjection::Project(origPnts, Npnts, copy);
  Float_t R, NR, y, z;
  for(Int_t i = 0; i<Npnts; i++) 
  {
    R = pnts[i].R();
    NR =R/(1+R*fDistortion);
    y = ( pnts[i].y > 0) ? NR : -NR;
    z = pnts[i].z;
    pnts[i].Set(z/(1+TMath::Abs(z*fDistortion)), y, 0);
  }
  return pnts;
}

//______________________________________________________________________________
Bool_t RhoZ::AcceptSegment(Vector& v1, Vector& v2) 
{
  return (v1.y*v2.y < 0) ? kFALSE : kTRUE;
}

//______________________________________________________________________________
Vector*  CircularFishEye::Project(Vector* origPnts, Int_t Npnts, Bool_t copy ) 
{
  Vector* pnts = NLTProjection::Project(origPnts, Npnts, copy);
  
  Float_t R, NR, phi;
  for(Int_t i = 0; i<Npnts; i++) {
    R = pnts[i].R();
    phi = pnts[i].Phi();
    NR = R/(1+R*fDistortion);
    pnts[i].Set(NR*TMath::Cos(phi), NR*TMath::Sin(phi), 0);
  }

  return pnts;
}

/**************************************************************************/
/**************************************************************************/
//______________________________________________________________________
// NLTProjector
//

ClassImp(NLTProjector)

NLTProjector::NLTProjector():
  TNamed("NLTProjector",""),
  fProjection(0),
  fEps(0.05),fIdxMap(0), 
  fNRPnts(0), fRPnts(0)
{
}

//______________________________________________________________________________
NLTProjector::~NLTProjector()
{
  CleanUp();

  if(fProjection) delete fProjection;
}

//______________________________________________________________________________
void NLTProjector::CleanUp()
{
  if(fIdxMap)
  {
    delete [] fIdxMap; 
    fIdxMap = 0;
  }
  if(fRPnts)
  {
    delete [] fRPnts; 
    fNRPnts  = 0;
    fRPnts = 0;
  }
}

/**************************************************************************/
void NLTProjector::SetProjection(NLTProjection::Type_e type, Float_t distort)
{
  static const Exc_t eH("NLTProjector::SetProjection ");
  switch (type)
  {
    case NLTProjection::RhoZ:
    {
      fProjection  = new RhoZ();
      break;
    }
    case NLTProjection::CFishEye:
    {
      fProjection  = new CircularFishEye();
      break;
    }
    default:
      throw(eH + "projection type not valid.");
      break;
  }
  fProjection->fDistortion = distort;
}

//______________________________________________________________________________
void NLTProjector::DumpBuffer(TBuffer3D* buff)
{
  Int_t* bpols = buff->fPols;
  
  for(UInt_t pi = 0; pi< buff->NbPols(); pi++) 
  {
    UInt_t Nseg = bpols[1]; 
    printf("%d polygon of %d has %d segments \n", pi,buff->NbPols(),Nseg);
    
    Int_t* seg =  &bpols[2];
    for(UInt_t a=0; a<Nseg; a++)
    {
      Int_t a1 = buff->fSegs[3*seg[a]+ 1];
      Int_t a2 = buff->fSegs[3*seg[a]+ 2];
      printf("(%d, %d) \n", a1, a2);
      printf("ORIG points :(%f, %f, %f)  (%f, %f, %f)\n", 
      	     buff->fPnts[3*a1],buff->fPnts[3*a1+1], buff->fPnts[3*a1+2],
      	     buff->fPnts[3*a2],buff->fPnts[3*a2+1], buff->fPnts[3*a2+2]);
    }
    printf("\n");
    bpols += (Nseg+2);
  }
}

//______________________________________________________________________________
void NLTProjector::CheckPoint(Int_t idx, Float_t* bbox)
{
  if(fRPnts[idx].x < bbox[0]) bbox[0] = fRPnts[idx].x;   
  if(fRPnts[idx].x > bbox[1]) bbox[1] = fRPnts[idx].x;

  if(fRPnts[idx].y < bbox[2]) bbox[2] = fRPnts[idx].y;   
  if(fRPnts[idx].y > bbox[3]) bbox[3] = fRPnts[idx].y;

  if(fRPnts[idx].z < bbox[4]) bbox[4] = fRPnts[idx].z;   
  if(fRPnts[idx].z > bbox[5]) bbox[5] = fRPnts[idx].z;
}

//______________________________________________________________________________
Bool_t NLTProjector::IsFirstIdxHead(Int_t s0, Int_t s1, TBuffer3D* buff)
{
  Int_t v0 = buff->fSegs[3*s0 + 1];
  Int_t v2 = buff->fSegs[3*s1 + 1];
  Int_t v3 = buff->fSegs[3*s1 + 2];
  if(v0 != v2 && v0 != v3 )
    return kTRUE;
  else 
    return kFALSE;
}

//______________________________________________________________________________
void  NLTProjector::ReducePoints(Vector* pnts, Int_t N)
{
  fIdxMap   = new Int_t[N];  
  Int_t* ra = new Int_t[N];  // list of reduced vertices
  fNRPnts = 0;
  
  for(UInt_t v = 0; v < (UInt_t)N; v++)
  {
    fIdxMap[v] = -1;
    for(Int_t k = 0; k<fNRPnts; k++) 
    {
      if(pnts[v].SquareDistance(pnts[ra[k]]) < fEps*fEps)
      {
	fIdxMap[v] = k; 
	break;
      }
    } 
    // have not find a point inside epsilon, add new point in scaled array
    if(fIdxMap[v] == -1)
    {
      fIdxMap[v] = fNRPnts;
      ra[fNRPnts] = v;
      fNRPnts++;
    }
    // printf("(%f, %f) vertex map %d -> %d \n", pnts[v*2], pnts[v*2 + 1], v, fIdxMap[v]);
  }
  
  // create an array of scaled points
  fRPnts = new Vector[fNRPnts];
  for(Int_t i = 0; i<fNRPnts; i++)
    fRPnts[i].Set(pnts[ra[i]].x,  pnts[ra[i]].y,  pnts[ra[i]].z);
  
  delete [] ra;  
  //  printf("reduced %d points of %d\n", fNRPnts, N);
}

//______________________________________________________________________________
void  NLTProjector::MakePolygonsFromBP(TBuffer3D* buff, std::list<NLTPolygon>& pols)
{
  //  printf("START NLTProjector::MakePolygonsFromBP\n");
  Int_t* bpols = buff->fPols;
  for(UInt_t pi = 0; pi< buff->NbPols(); pi++) 
  {
    std::list<Int_t>  pp; // points in current polygon 
    UInt_t Nseg = bpols[1]; 
    Int_t* seg =  &bpols[2];

    // start idx in the fist segment depends of second segment 
    Int_t  tail, head;
    Bool_t h = IsFirstIdxHead(seg[0], seg[1], buff);
    if(h) {
      head = fIdxMap[buff->fSegs[3*seg[0] + 1]];
      tail = fIdxMap[buff->fSegs[3*seg[0] + 2]];
    }
    else {
      head = fIdxMap[buff->fSegs[3*seg[0] + 2]];
      tail = fIdxMap[buff->fSegs[3*seg[0] + 1]];
    }
    Float_t bbox[] = {0., 0., 0., 0., 0., 0.};
    pp.push_back(head);
    CheckPoint(head, bbox);
    // printf("start idx head %d, tail %d\n", head, tail);
    
    std::list<Seg> segs;  
    for(UInt_t s=1; s< Nseg; s++)
      segs.push_back(Seg(buff->fSegs[3*seg[s] + 1],buff->fSegs[3*seg[s] + 2]));
    Bool_t accepted = kFALSE; 
    for(std::list<Seg>::iterator it = segs.begin(); it != segs.end(); it++ )
    { 
      Int_t mv1 = fIdxMap[(*it).v1];
      Int_t mv2 = fIdxMap[(*it).v2];         
      accepted = fProjection->AcceptSegment(fRPnts[mv1], fRPnts[mv2]);
      if(accepted == kFALSE)
      {
	pp.clear();
	break;
      }	  
      if(tail != pp.back()) 
      {
	pp.push_back(tail);
	CheckPoint(tail, bbox);
      }
      tail = (mv1 == tail) ? mv2 :mv1;
    }

    // render class loops indices last and first should not be equal
    if(accepted && pp.front() == pp.back())
      pp.pop_front();
    

    if( pp.size()>2 && (bbox[1]-bbox[0])>fEps && (bbox[3]-bbox[2])>fEps)
    {
      Bool_t duplicate = kFALSE;
      for (std::list<NLTPolygon>::iterator poi = pols.begin(); poi!= pols.end(); poi++)
      {
        NLTPolygon P = *poi;
        if(pp.size() != (UInt_t)P.fNPnts) 
	  continue;      
	Int_t matched = 0;
	for (std::list<Int_t>::iterator u = pp.begin(); u!= pp.end(); u++) {
          if (*u == P.fPnts[*u]) matched++;
          else continue;
        }
        if (matched == P.fNPnts) {
          duplicate = kTRUE;
          break;
        }
      }

      if(duplicate == kFALSE) 
      {
	Int_t* pv = new Int_t[pp.size()];
	Int_t count=0;
	// printf("%d NLTPolygon points %d \n", pols.size(), pp.size());
        // printf("bbox (%f, %f) (%f, %f)\n", bbox[1], bbox[0], bbox[3], bbox[2]);
	for( std::list<Int_t>::iterator u = pp.begin(); u!= pp.end(); u++){
	  pv[count] = *u;
	  count++;
	}
	pols.push_back(NLTPolygon(pp.size(), pv));
      }
    } // end of NLTPolygon
    bpols += (Nseg+2);
  }
}// MakePolygonsFromBP

//______________________________________________________________________________
void  NLTProjector::MakePolygonsFromBS(TBuffer3D* buff, std::list<NLTPolygon>& pols)
{
  // create your own list of segments according to reduced and projected points
  std::list<Seg> segs;  
  std::list<Seg>::iterator it;
  for(UInt_t s=0; s< buff->NbSegs(); s++)
  {
    Bool_t duplicate = kFALSE;
    Int_t vo1, vo2;   // idx from buff segment
    Int_t vor1, vor2; // mapped idx 
    vo1 =  buff->fSegs[3*s + 1];
    vo2 =  buff->fSegs[3*s + 2]; //... skip color info
    vor1 = fIdxMap[vo1];
    vor2 = fIdxMap[vo2];
    if(vor1 == vor2) continue;
    // check duplicate
    for(it = segs.begin(); it != segs.end(); it++ ){
      Int_t vv1 = (*it).v1;
      Int_t vv2 = (*it).v2;
      if((vv1 == vor1 && vv2 == vor2 )||(vv1 == vor2 && vv2 == vor1 )){
	duplicate = kTRUE;
	continue;
      }
    }
    if(duplicate == kFALSE && fProjection->AcceptSegment(fRPnts[vor1], fRPnts[vor2]))
    {
      segs.push_back(Seg(vor1, vor2));
    }
  }

  while(segs.empty() == kFALSE)
  {
    // printf("Start building polygon %d from %d segments in POOL \n", pols.size(), segs.size());
    std::list<Int_t> pp; // points in current polygon
    pp.push_back(segs.front().v1);
    Int_t tail = segs.front().v2;
    segs.pop_front();
    Bool_t match = kTRUE;
    while(match && segs.empty() == kFALSE)
    {
      // printf("second loop search tail %d \n",tail);
      for(It_t k=segs.begin(); k!=segs.end(); ++k){
	Int_t cv1 = (*k).v1;
	Int_t cv2 = (*k).v2;
	if( cv1 == tail || cv2 == tail){
	  // printf("found point %d  in %d,%d  \n", tail, cv1, cv2);
	  pp.push_back(tail);
	  if(cv1 == tail) 
	    tail = cv2;
	  else 
	    tail = cv1;

	  It_t to_erase = k--;
	  segs.erase(to_erase);
	  match = kTRUE;
	  break;
	}
	else 
	{
	  match = kFALSE;
	}
      } // end for loop in the segment pool
      if(tail ==  pp.front())
	break;
    };

    if(pp.size() > 2){      
      // printf("CLOSING polygon %d with %d points \n", pols.size(),pp.size());
      Int_t* pv = new Int_t[pp.size()];
      Int_t count=0;
      for( std::list<Int_t>::iterator u = pp.begin(); u!= pp.end(); u++){
	pv[count] = *u;
	count++;
      }
      pols.push_back(NLTPolygon(pp.size(), pv));
    }
  } // while segment pool empty
}//MakePolygonsFromBS

 //______________________________________________________________________________
NLTPolygonSet*  NLTProjector::ProjectGeoShape(TBuffer3D* buff, Int_t useBuffPols)
{ 
  // DumpBuffer(buff);  
  // rewrite from Double_t to Vector array
  Vector*  pnts  = new Vector[buff->NbPnts()];
  for(Int_t i = 0; i<(Int_t)buff->NbPnts(); i++) 
    pnts[i].Set(buff->fPnts[3*i],buff->fPnts[3*i+1], buff->fPnts[3*i+2]);

  fProjection->Project(pnts, buff->NbPnts(), kFALSE);
  ReducePoints(pnts, buff->NbPnts());

  // build polygons
  std::list<NLTPolygon>   pols;
  if(useBuffPols == -1) 
  { 
    MakePolygonsFromBP(buff, pols); 
    if(pols.empty())
    {
      // printf("useBuffPols == -1 call MakePolygonsFromBS \n");
      MakePolygonsFromBS(buff, pols);
    }
  }
  if(useBuffPols == 1)
  {
    MakePolygonsFromBP(buff, pols); 
  }
  if(useBuffPols == 0) 
  {
    MakePolygonsFromBS(buff, pols); 
  }
 
  NLTPolygonSet* ps = new NLTPolygonSet();
  if(pols.empty() == kFALSE) 
  {
    // points
    ps->SetPoints(fRPnts,fNRPnts);
    // polygons
    NLTPolygon* nltp = new NLTPolygon[pols.size()];
    Int_t pc = 0;
    for( std::list<NLTPolygon>::iterator pi = pols.begin(); pi!= pols.end(); pi++)
    {
      nltp[pc].fNPnts = (*pi).fNPnts;
      nltp[pc].fPnts = (*pi).fPnts;
      pc++;
    }
    ps->SetPolygons(nltp, pols.size());
    // ps->Dump();
  }
  fRPnts = 0;
  CleanUp();
  return ps;
}

//______________________________________________________________________________
void NLTProjector::ProjectPointSet(PointSet* orig)
{
  // keep original rewrite points into Vector array
  Float_t* op = orig->GetP();
  Vector* vp = new Vector[orig->GetN()];
  for( Int_t k=0; k<orig->GetN(); k++)
    vp[k].Set(op[3*k], op[3*k+1], op[3*k+2]);
  Vector* pnts = fProjection->Project(vp, orig->GetN(), kFALSE);
  ReducePoints(pnts, orig->GetN());
  orig->Reset(fNRPnts);
  for(Int_t i = 0; i< fNRPnts; i++){
    orig->SetPoint(i, fRPnts[i].x, fRPnts[i].y, fRPnts[3].z);
  }

  CleanUp();
}
  
//______________________________________________________________________________
void  NLTProjector::ProjectTrackList(TrackList* orig_tl)
{
  // fill a track list from original track list 
  Reve::RenderElement::List_i i = orig_tl->BeginChildren();
  Int_t Nc =  orig_tl->GetNChildren();
  for(Int_t c=0; c<Nc; c++)
  {
    // printf("project track with %d points \n", orig->GetN());
    Track* orig = dynamic_cast<Track*>(*i);
    Int_t start = 0;
    // keep original copy points to Vector array
    Float_t* op = orig->GetP();
    Vector* pnts = new Vector[orig->GetN()];
    for( Int_t k=0; k<orig->GetN(); k++)
      pnts[k].Set(op[3*k], op[3*k+1], op[3*k+2]);

    for( Int_t k=0; k<orig->GetN()-1; k++)
    {
      if(fProjection->AcceptSegment(pnts[k], pnts[(k+1)]) == kFALSE){
        //printf("split track at %d\n", k);
	Vector* tp = fProjection->Project(&pnts[start], k-start, kFALSE);
        ReducePoints(tp, k-start);
	Track* nt = new Track();  
        nt->Reset(fNRPnts);
	nt->SetMainColor(orig->GetMainColor());
        for(Int_t j = 0; j<fNRPnts; j++)
          nt->SetPoint(j, fRPnts[j].x, fRPnts[j].y, fRPnts[j].z);
        orig_tl->AddElement(nt);
	CleanUp();
        start = k+1;
      }
    }
    if(start != orig->GetN()-1) {
      // printf("finish last track %d %d\n", start,orig->GetN()-1);
      Track* nt = new Track();
      Vector* tp = fProjection->Project(&pnts[start], orig->GetN()-start, kFALSE);
      ReducePoints(tp, orig->GetN()-start);
      nt->Reset(fNRPnts);
      nt->SetMainColor(orig->GetMainColor());
      for(Int_t j = 0; j<fNRPnts; j++)
	nt->SetPoint(j, fRPnts[j].x, fRPnts[j].y, fRPnts[j].z);
      orig_tl->AddElement(nt);
      CleanUp();
    }
    // destroy unprojected
    ++i;
  } // end loop tracks

  // remove original tracks
  for(Int_t c=0; c<Nc; c++)
  {
    i = orig_tl->BeginChildren();
    (*i)->Destroy();
  }
}
Commit	Line	Data
debf9f47	1	#include "NLTProjector.h"
	2	#include "RGTopFrame.h"
	3	#include "NLTPolygonSet.h"
	4	#include "PODs.h"
	5	#include "PointSet.h"
	6	#include "Track.h"
	7
	8	#include "TBuffer3D.h"
	9	#include "TColor.h"
	10	#include "TPointSet3D.h"
	11	#include "TGeoManager.h"
	12	#include "TGeoMatrix.h"
	13
	14	#include <list>
	15
	16	using namespace Reve;
	17
	18	namespace {
	19	struct Seg {
	20	Int_t v1;
	21	Int_t v2;
	22
	23	Seg(Int_t i1=-1, Int_t i2=-1):v1(i1), v2(i2){};
	24	};
	25
	26	typedef std::list<Seg>::iterator It_t;
	27	}
	28
	29	//______________________________________________________________________________
	30	Vector* NLTProjection::Project(Vector* origPnts, Int_t Npnts, Bool_t copy)
	31	{
	32	Vector* pnts = 0;
	33	if(copy)
	34	{
	35	pnts = (Vector* )malloc(Npnts*sizeof(Vector));
	36	memcpy(pnts, origPnts, Npnts*sizeof(Vector));
	37	}
	38	else
	39	{
	40	pnts = origPnts;
	41	}
	42	return pnts;
	43	}
	44
	45	//______________________________________________________________________________
	46	Vector* RhoZ::Project(Vector* origPnts, Int_t Npnts, Bool_t copy)
	47	{
	48	Vector* pnts = NLTProjection::Project(origPnts, Npnts, copy);
	49	Float_t R, NR, y, z;
	50	for(Int_t i = 0; i<Npnts; i++)
	51	{
	52	R = pnts[i].R();
	53	NR =R/(1+R*fDistortion);
	54	y = ( pnts[i].y > 0) ? NR : -NR;
	55	z = pnts[i].z;
	56	pnts[i].Set(z/(1+TMath::Abs(z*fDistortion)), y, 0);
	57	}
	58	return pnts;
	59	}
	60
	61	//______________________________________________________________________________
	62	Bool_t RhoZ::AcceptSegment(Vector& v1, Vector& v2)
	63	{
	64	return (v1.y*v2.y < 0) ? kFALSE : kTRUE;
65	}
66
67	//______________________________________________________________________________
68	Vector* CircularFishEye::Project(Vector* origPnts, Int_t Npnts, Bool_t copy )
69	{
70	Vector* pnts = NLTProjection::Project(origPnts, Npnts, copy);
71
72	Float_t R, NR, phi;
73	for(Int_t i = 0; i<Npnts; i++) {
74	R = pnts[i].R();
75	phi = pnts[i].Phi();
76	NR = R/(1+R*fDistortion);
77	pnts[i].Set(NRTMath::Cos(phi), NRTMath::Sin(phi), 0);
78	}
79
80	return pnts;
81	}
82
83	/**************************************************************************/
84	/**************************************************************************/
85	//______________________________________________________________________
86	// NLTProjector
87	//
88
89	ClassImp(NLTProjector)
90
91	NLTProjector::NLTProjector():
92	TNamed("NLTProjector",""),
93	fProjection(0),
94	fEps(0.05),fIdxMap(0),
95	fNRPnts(0), fRPnts(0)
96	{
97	}
98
99	//______________________________________________________________________________
100	NLTProjector::~NLTProjector()
101	{
102	CleanUp();
103
104	if(fProjection) delete fProjection;
105	}
106
107	//______________________________________________________________________________
108	void NLTProjector::CleanUp()
109	{
110	if(fIdxMap)
111	{
112	delete [] fIdxMap;
113	fIdxMap = 0;
114	}
115	if(fRPnts)
116	{
117	delete [] fRPnts;
118	fNRPnts = 0;
119	fRPnts = 0;
120	}
121	}
122
123	/**************************************************************************/
124	void NLTProjector::SetProjection(NLTProjection::Type_e type, Float_t distort)
125	{
126	static const Exc_t eH("NLTProjector::SetProjection ");
127	switch (type)
128	{
129	case NLTProjection::RhoZ:
130	{
131	fProjection = new RhoZ();
132	break;
133	}
134	case NLTProjection::CFishEye:
135	{
136	fProjection = new CircularFishEye();
137	break;
138	}
139	default:
140	throw(eH + "projection type not valid.");
141	break;
142	}
143	fProjection->fDistortion = distort;
144	}
145
146	//______________________________________________________________________________
147	void NLTProjector::DumpBuffer(TBuffer3D* buff)
148	{
149	Int_t* bpols = buff->fPols;
150
151	for(UInt_t pi = 0; pi< buff->NbPols(); pi++)
152	{
153	UInt_t Nseg = bpols[1];
154	printf("%d polygon of %d has %d segments \n", pi,buff->NbPols(),Nseg);
155
156	Int_t* seg = &bpols[2];
157	for(UInt_t a=0; a<Nseg; a++)
158	{
159	Int_t a1 = buff->fSegs[3*seg[a]+ 1];
160	Int_t a2 = buff->fSegs[3*seg[a]+ 2];
161	printf("(%d, %d) \n", a1, a2);
162	printf("ORIG points :(%f, %f, %f) (%f, %f, %f)\n",
163	buff->fPnts[3a1],buff->fPnts[3a1+1], buff->fPnts[3*a1+2],
164	buff->fPnts[3a2],buff->fPnts[3a2+1], buff->fPnts[3*a2+2]);
165	}
166	printf("\n");
167	bpols += (Nseg+2);
168	}
169	}
170
171	//______________________________________________________________________________
172	void NLTProjector::CheckPoint(Int_t idx, Float_t* bbox)
173	{
174	if(fRPnts[idx].x < bbox[0]) bbox[0] = fRPnts[idx].x;
175	if(fRPnts[idx].x > bbox[1]) bbox[1] = fRPnts[idx].x;
176
177	if(fRPnts[idx].y < bbox[2]) bbox[2] = fRPnts[idx].y;
178	if(fRPnts[idx].y > bbox[3]) bbox[3] = fRPnts[idx].y;
179
180	if(fRPnts[idx].z < bbox[4]) bbox[4] = fRPnts[idx].z;
181	if(fRPnts[idx].z > bbox[5]) bbox[5] = fRPnts[idx].z;
182	}
183
184	//______________________________________________________________________________
185	Bool_t NLTProjector::IsFirstIdxHead(Int_t s0, Int_t s1, TBuffer3D* buff)
186	{
187	Int_t v0 = buff->fSegs[3*s0 + 1];
188	Int_t v2 = buff->fSegs[3*s1 + 1];
189	Int_t v3 = buff->fSegs[3*s1 + 2];
190	if(v0 != v2 && v0 != v3 )
191	return kTRUE;
192	else
193	return kFALSE;
194	}
195
196	//______________________________________________________________________________
197	void NLTProjector::ReducePoints(Vector* pnts, Int_t N)
198	{
199	fIdxMap = new Int_t[N];
200	Int_t* ra = new Int_t[N]; // list of reduced vertices
201	fNRPnts = 0;
202
203	for(UInt_t v = 0; v < (UInt_t)N; v++)
204	{
205	fIdxMap[v] = -1;
206	for(Int_t k = 0; k<fNRPnts; k++)
207	{
208	if(pnts[v].SquareDistance(pnts[ra[k]]) < fEps*fEps)
209	{
210	fIdxMap[v] = k;
211	break;
212	}
213	}
214	// have not find a point inside epsilon, add new point in scaled array
215	if(fIdxMap[v] == -1)
216	{
217	fIdxMap[v] = fNRPnts;
218	ra[fNRPnts] = v;
219	fNRPnts++;
220	}
221	// printf("(%f, %f) vertex map %d -> %d \n", pnts[v2], pnts[v2 + 1], v, fIdxMap[v]);
222	}
223
224	// create an array of scaled points
225	fRPnts = new Vector[fNRPnts];
226	for(Int_t i = 0; i<fNRPnts; i++)
227	fRPnts[i].Set(pnts[ra[i]].x, pnts[ra[i]].y, pnts[ra[i]].z);
228
229	delete [] ra;
230	// printf("reduced %d points of %d\n", fNRPnts, N);
231	}
232
233	//______________________________________________________________________________
234	void NLTProjector::MakePolygonsFromBP(TBuffer3D* buff, std::list<NLTPolygon>& pols)
235	{
236	// printf("START NLTProjector::MakePolygonsFromBP\n");
237	Int_t* bpols = buff->fPols;
238	for(UInt_t pi = 0; pi< buff->NbPols(); pi++)
239	{
240	std::list<Int_t> pp; // points in current polygon
241	UInt_t Nseg = bpols[1];
242	Int_t* seg = &bpols[2];
243
244	// start idx in the fist segment depends of second segment
245	Int_t tail, head;
246	Bool_t h = IsFirstIdxHead(seg[0], seg[1], buff);
247	if(h) {
248	head = fIdxMap[buff->fSegs[3*seg[0] + 1]];
249	tail = fIdxMap[buff->fSegs[3*seg[0] + 2]];
250	}
251	else {
252	head = fIdxMap[buff->fSegs[3*seg[0] + 2]];
253	tail = fIdxMap[buff->fSegs[3*seg[0] + 1]];
254	}
255	Float_t bbox[] = {0., 0., 0., 0., 0., 0.};
256	pp.push_back(head);
257	CheckPoint(head, bbox);
258	// printf("start idx head %d, tail %d\n", head, tail);
259
260	std::list<Seg> segs;
261	for(UInt_t s=1; s< Nseg; s++)
262	segs.push_back(Seg(buff->fSegs[3seg[s] + 1],buff->fSegs[3seg[s] + 2]));
263	Bool_t accepted = kFALSE;
264	for(std::list<Seg>::iterator it = segs.begin(); it != segs.end(); it++ )
265	{
266	Int_t mv1 = fIdxMap[(*it).v1];
267	Int_t mv2 = fIdxMap[(*it).v2];
268	accepted = fProjection->AcceptSegment(fRPnts[mv1], fRPnts[mv2]);
269	if(accepted == kFALSE)
270	{
271	pp.clear();
272	break;
273	}
274	if(tail != pp.back())
275	{
276	pp.push_back(tail);
277	CheckPoint(tail, bbox);
278	}
279	tail = (mv1 == tail) ? mv2 :mv1;
280	}
281
282	// render class loops indices last and first should not be equal
283	if(accepted && pp.front() == pp.back())
284	pp.pop_front();
285
286
287	if( pp.size()>2 && (bbox[1]-bbox[0])>fEps && (bbox[3]-bbox[2])>fEps)
288	{
289	Bool_t duplicate = kFALSE;
f36b73bd	290	for (std::list<NLTPolygon>::iterator poi = pols.begin(); poi!= pols.end(); poi++)
debf9f47	291	{
f36b73bd	292	NLTPolygon P = *poi;
debf9f47	293	if(pp.size() != (UInt_t)P.fNPnts)
	294	continue;
	295	Int_t matched = 0;
f36b73bd	296	for (std::list<Int_t>::iterator u = pp.begin(); u!= pp.end(); u++) {
f36b73bd	297	if (u == P.fPnts[u]) matched++;
debf9f47	298	else continue;
debf9f47	299	}
f36b73bd	300	if (matched == P.fNPnts) {
debf9f47	301	duplicate = kTRUE;
	302	break;
	303	}
	304	}
	305
	306	if(duplicate == kFALSE)
	307	{
	308	Int_t* pv = new Int_t[pp.size()];
	309	Int_t count=0;
	310	// printf("%d NLTPolygon points %d \n", pols.size(), pp.size());
	311	// printf("bbox (%f, %f) (%f, %f)\n", bbox[1], bbox[0], bbox[3], bbox[2]);
	312	for( std::list<Int_t>::iterator u = pp.begin(); u!= pp.end(); u++){
	313	pv[count] = *u;
	314	count++;
	315	}
	316	pols.push_back(NLTPolygon(pp.size(), pv));
	317	}
	318	} // end of NLTPolygon
	319	bpols += (Nseg+2);
	320	}
	321	}// MakePolygonsFromBP
	322
	323	//______________________________________________________________________________
	324	void NLTProjector::MakePolygonsFromBS(TBuffer3D* buff, std::list<NLTPolygon>& pols)
	325	{
	326	// create your own list of segments according to reduced and projected points
	327	std::list<Seg> segs;
	328	std::list<Seg>::iterator it;
	329	for(UInt_t s=0; s< buff->NbSegs(); s++)
	330	{
	331	Bool_t duplicate = kFALSE;
	332	Int_t vo1, vo2; // idx from buff segment
	333	Int_t vor1, vor2; // mapped idx
	334	vo1 = buff->fSegs[3*s + 1];
	335	vo2 = buff->fSegs[3*s + 2]; //... skip color info
	336	vor1 = fIdxMap[vo1];
	337	vor2 = fIdxMap[vo2];
	338	if(vor1 == vor2) continue;
	339	// check duplicate
	340	for(it = segs.begin(); it != segs.end(); it++ ){
	341	Int_t vv1 = (*it).v1;
	342	Int_t vv2 = (*it).v2;
	343	if((vv1 == vor1 && vv2 == vor2 )\|\|(vv1 == vor2 && vv2 == vor1 )){
	344	duplicate = kTRUE;
	345	continue;
	346	}
	347	}
	348	if(duplicate == kFALSE && fProjection->AcceptSegment(fRPnts[vor1], fRPnts[vor2]))
	349	{
	350	segs.push_back(Seg(vor1, vor2));
	351	}
	352	}
	353
	354	while(segs.empty() == kFALSE)
	355	{
	356	// printf("Start building polygon %d from %d segments in POOL \n", pols.size(), segs.size());
	357	std::list<Int_t> pp; // points in current polygon
	358	pp.push_back(segs.front().v1);
	359	Int_t tail = segs.front().v2;
	360	segs.pop_front();
	361	Bool_t match = kTRUE;
	362	while(match && segs.empty() == kFALSE)
	363	{
	364	// printf("second loop search tail %d \n",tail);
365	for(It_t k=segs.begin(); k!=segs.end(); ++k){
366	Int_t cv1 = (*k).v1;
367	Int_t cv2 = (*k).v2;
368	if( cv1 == tail \|\| cv2 == tail){
369	// printf("found point %d in %d,%d \n", tail, cv1, cv2);
370	pp.push_back(tail);
371	if(cv1 == tail)
372	tail = cv2;
373	else
374	tail = cv1;
375
376	It_t to_erase = k--;
377	segs.erase(to_erase);
378	match = kTRUE;
379	break;
380	}
381	else
382	{
383	match = kFALSE;
384	}
385	} // end for loop in the segment pool
386	if(tail == pp.front())
387	break;
388	};
389
390	if(pp.size() > 2){
391	// printf("CLOSING polygon %d with %d points \n", pols.size(),pp.size());
392	Int_t* pv = new Int_t[pp.size()];
393	Int_t count=0;
394	for( std::list<Int_t>::iterator u = pp.begin(); u!= pp.end(); u++){
395	pv[count] = *u;
396	count++;
397	}
398	pols.push_back(NLTPolygon(pp.size(), pv));
399	}
400	} // while segment pool empty
401	}//MakePolygonsFromBS
402
403	//______________________________________________________________________________
404	NLTPolygonSet* NLTProjector::ProjectGeoShape(TBuffer3D* buff, Int_t useBuffPols)
405	{
406	// DumpBuffer(buff);
407	// rewrite from Double_t to Vector array
408	Vector* pnts = new Vector[buff->NbPnts()];
409	for(Int_t i = 0; i<(Int_t)buff->NbPnts(); i++)
410	pnts[i].Set(buff->fPnts[3i],buff->fPnts[3i+1], buff->fPnts[3*i+2]);
411
412	fProjection->Project(pnts, buff->NbPnts(), kFALSE);
413	ReducePoints(pnts, buff->NbPnts());
414
415	// build polygons
416	std::list<NLTPolygon> pols;
417	if(useBuffPols == -1)
418	{
419	MakePolygonsFromBP(buff, pols);
420	if(pols.empty())
421	{
422	// printf("useBuffPols == -1 call MakePolygonsFromBS \n");
423	MakePolygonsFromBS(buff, pols);
424	}
425	}
426	if(useBuffPols == 1)
427	{
428	MakePolygonsFromBP(buff, pols);
429	}
430	if(useBuffPols == 0)
431	{
432	MakePolygonsFromBS(buff, pols);
433	}
434
435	NLTPolygonSet* ps = new NLTPolygonSet();
436	if(pols.empty() == kFALSE)
437	{
438	// points
439	ps->SetPoints(fRPnts,fNRPnts);
440	// polygons
441	NLTPolygon* nltp = new NLTPolygon[pols.size()];
442	Int_t pc = 0;
443	for( std::list<NLTPolygon>::iterator pi = pols.begin(); pi!= pols.end(); pi++)
444	{
445	nltp[pc].fNPnts = (*pi).fNPnts;
446	nltp[pc].fPnts = (*pi).fPnts;
447	pc++;
448	}
449	ps->SetPolygons(nltp, pols.size());
450	// ps->Dump();
451	}
452	fRPnts = 0;
453	CleanUp();
454	return ps;
455	}
456
457	//______________________________________________________________________________
458	void NLTProjector::ProjectPointSet(PointSet* orig)
459	{
460	// keep original rewrite points into Vector array
461	Float_t* op = orig->GetP();
462	Vector* vp = new Vector[orig->GetN()];
463	for( Int_t k=0; k<orig->GetN(); k++)
464	vp[k].Set(op[3k], op[3k+1], op[3*k+2]);
465	Vector* pnts = fProjection->Project(vp, orig->GetN(), kFALSE);
466	ReducePoints(pnts, orig->GetN());
467	orig->Reset(fNRPnts);
468	for(Int_t i = 0; i< fNRPnts; i++){
469	orig->SetPoint(i, fRPnts[i].x, fRPnts[i].y, fRPnts[3].z);
470	}
471
472	CleanUp();
473	}
474
475	//______________________________________________________________________________
476	void NLTProjector::ProjectTrackList(TrackList* orig_tl)
477	{
478	// fill a track list from original track list
479	Reve::RenderElement::List_i i = orig_tl->BeginChildren();
480	Int_t Nc = orig_tl->GetNChildren();
481	for(Int_t c=0; c<Nc; c++)
482	{
483	// printf("project track with %d points \n", orig->GetN());
484	Track* orig = dynamic_cast<Track>(i);
485	Int_t start = 0;
486	// keep original copy points to Vector array
487	Float_t* op = orig->GetP();
488	Vector* pnts = new Vector[orig->GetN()];
489	for( Int_t k=0; k<orig->GetN(); k++)
490	pnts[k].Set(op[3k], op[3k+1], op[3*k+2]);
491
492	for( Int_t k=0; k<orig->GetN()-1; k++)
493	{
494	if(fProjection->AcceptSegment(pnts[k], pnts[(k+1)]) == kFALSE){
495	//printf("split track at %d\n", k);
496	Vector* tp = fProjection->Project(&pnts[start], k-start, kFALSE);
497	ReducePoints(tp, k-start);
498	Track* nt = new Track();
499	nt->Reset(fNRPnts);
500	nt->SetMainColor(orig->GetMainColor());
501	for(Int_t j = 0; j<fNRPnts; j++)
502	nt->SetPoint(j, fRPnts[j].x, fRPnts[j].y, fRPnts[j].z);
503	orig_tl->AddElement(nt);
504	CleanUp();
505	start = k+1;
506	}
507	}
508	if(start != orig->GetN()-1) {
509	// printf("finish last track %d %d\n", start,orig->GetN()-1);
510	Track* nt = new Track();
511	Vector* tp = fProjection->Project(&pnts[start], orig->GetN()-start, kFALSE);
512	ReducePoints(tp, orig->GetN()-start);
513	nt->Reset(fNRPnts);
514	nt->SetMainColor(orig->GetMainColor());
515	for(Int_t j = 0; j<fNRPnts; j++)
516	nt->SetPoint(j, fRPnts[j].x, fRPnts[j].y, fRPnts[j].z);
517	orig_tl->AddElement(nt);
518	CleanUp();
519	}
520	// destroy unprojected
521	++i;
522	} // end loop tracks
523
524	// remove original tracks
525	for(Int_t c=0; c<Nc; c++)
526	{
527	i = orig_tl->BeginChildren();
528	(*i)->Destroy();
529	}
530	}