EVE/Reve/NLTProjector.cxx

   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(NR*TMath::Cos(phi), NR*TMath::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[3*a1],buff->fPnts[3*a1+1], buff->fPnts[3*a1+2],
 164              buff->fPnts[3*a2],buff->fPnts[3*a2+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[v*2], pnts[v*2 + 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[3*seg[s] + 1],buff->fSegs[3*seg[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;
 290       for (std::list<NLTPolygon>::iterator poi = pols.begin(); poi!= pols.end(); poi++)
 291       {
 292         NLTPolygon P = *poi;
 293         if(pp.size() != (UInt_t)P.fNPnts)
 294           continue;
 295         Int_t matched = 0;
 296         for (std::list<Int_t>::iterator u = pp.begin(); u!= pp.end(); u++) {
 297           if (*u == P.fPnts[*u]) matched++;
 298           else continue;
 299         }
 300         if (matched == P.fNPnts) {
 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[3*i],buff->fPnts[3*i+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[3*k], op[3*k+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[3*k], op[3*k+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 }