EdbViewMatch Class Reference

#include <EdbViewMatch.h>

Inheritance diagram for EdbViewMatch:

Collaboration diagram for EdbViewMatch:

Public Member Functions
void	AddCluster (float x, float y, float z, float xv, float yv, int view, int frame)
void	CalculateCorr ()
void	CalculateGrRef ()
void	CutGrRef ()
void	DrawCorrMap ()
TNtuple *	DumpGr (const char *name)
	EdbViewMatch ()
void	GenerateCorrectionMatrix (const char *addfile)
void	InitCorrMap ()
void	InitGMap ()
void	MakeDistortionMap (const char *fname, TEnv &env, const char *addfile=0)
void	Print ()
void	SetPar (TEnv &env)
void	SetPixelSize (float xpix, float ypix)
virtual	~EdbViewMatch ()

Public Attributes
float	eAreaMax
float	eAreaMin
bool	eDumpGr
TFile *	eOutputFile
float	eVolumeMax
float	eVolumeMin

Private Attributes
TClonesArray	eCl
float	eCorrectionMatrixStepX
float	eCorrectionMatrixStepY
EdbCell2	eCorrMap
EdbCell2	eGMap
TClonesArray	eGr
int	eNClMin
int	eNXpix
int	eNYpix
float	eR2CenterMax
float	eRmax
float	eXpix
float	eYpix

Constructor & Destructor Documentation

EdbViewMatch()

EdbViewMatch::EdbViewMatch

(

)

~EdbViewMatch()

virtual EdbViewMatch::~EdbViewMatch ( )

inlinevirtual

{}

Member Function Documentation

AddCluster()

void EdbViewMatch::AddCluster	(	float	x,
		float	y,
		float	z,
		float	xv,
		float	yv,
		int	view,
		int	frame
	)

{
  int ind = eCl.GetEntriesFast();
  EdbClMatch *c = new(eCl[ind]) EdbClMatch(x,y,z, xv,yv, view,frame );
 
  float v[2]={ x+xv,  y+yv };
  int  ir[2]={1,1};
  TObjArray arr;
  EdbSegment *s0 = 0;
  int n = eGMap.SelectObjectsC( v, ir, arr);
  if(n) { 
   float r, rmin=2.*eRmax;
    for(int i=0; i<n; i++) {
      EdbSegment *s = (EdbSegment *)arr.At(i);
      r = Sqrt( (s->GetX0()-v[0])*(s->GetX0()-v[0]) + (s->GetY0()-v[1])*(s->GetY0()-v[1]) );
      if(r<eRmax) if(r<rmin) { rmin=r; s0 = s; }
    }
  }
  if(s0) { s0->AddElement(c); s0->SetX0(v[0]); s0->SetY0(v[1]); s0->SetZ0(z); }
  else  {
    int inds = eGr.GetEntriesFast();
    EdbSegment *s = new(eGr[inds]) EdbSegment(v[0],v[1],z,0,0);
    s->AddElement(c);
    eGMap.AddObject(v, s);
  }
}

CalculateCorr()

void EdbViewMatch::CalculateCorr

(

)

{
  int ngr = eGr.GetEntries();
  Log(2,"EdbViewMatch::CalculateCorr","%d grains used",ngr);
  for(int i=0; i<ngr; i++) {
    EdbSegment *s = ((EdbSegment*)eGr.UncheckedAt(i));
    int nc = s->GetNelements();                                  if(!nc) continue;
    float x0 = s->GetX0(), y0= s->GetY0();
    for( int ic=0; ic<nc; ic++ ) {
      EdbClMatch *c = (EdbClMatch *)(s->GetElements()->UncheckedAt(ic));
      float dx = c->eX+c->eXv - x0;
      float dy = c->eY+c->eYv - y0;
      TArrayD *dxy = (TArrayD*)(eCorrMap.GetObject(c->eX, c->eY, 0));
      if(!dxy) { Log(1,"EdbViewMatch::CalculateCorr","WARNING: null dxy:  ic=%d i=%d   cx,cy: %f %f",ic,i, c->eX, c->eY);       continue; }
      (*dxy)[0] += dx;
      (*dxy)[1] += dy;
      eCorrMap.Fill(c->eX, c->eY);
    }
  }
  
  int nc=eCorrMap.Ncell();
  Log(2,"EdbViewMatch::CalculateCorr","%d cells in the map",nc);
  for( int i=0; i<nc; i++ ) {
    TArrayD *dxy = (TArrayD*)( eCorrMap.GetObject(i, 0 ) );    if(!dxy) continue;
    int w = eCorrMap.Bin(i);                                   if(!w) continue;
    (*dxy)[0] /= w;
    (*dxy)[1] /= w;
  }
 
}

CalculateGrRef()

void EdbViewMatch::CalculateGrRef

(

)

{
  // take as grain reference position of the most central cluster
  
  int ngr = eGr.GetEntries();
  for(int i=0; i<ngr; i++) {
    EdbSegment *s = ((EdbSegment*)eGr.UncheckedAt(i));
    int nc = s->GetNelements();
    //if(nc<eNClMin)        { s->GetElements()->Clear(); continue; }
    float rmin=10000., r;
    for( int ic=0; ic<nc; ic++ ) {
      EdbClMatch *c = (EdbClMatch *)(s->GetElements()->UncheckedAt(ic));
      r = Sqrt( c->eX*c->eX + c->eY*c->eY );                                     // distance to view center
      if(r<rmin) { rmin=r; s->SetX0(c->eXv+c->eX); s->SetY0(c->eYv+c->eY); s->SetDz(r); }
    }
    //if(rmin>eR2CenterMax) s->GetElements()->Clear();
 }
}

CutGrRef()

void EdbViewMatch::CutGrRef

(

)

{
  // remove marginal grains from calculation
  int ngr = eGr.GetEntries();
  for(int i=0; i<ngr; i++) {
    EdbSegment *s = ((EdbSegment*)eGr.UncheckedAt(i));
    int nc = s->GetNelements();
    if     (nc<eNClMin)              s->GetElements()->Clear();
    else if(s->GetDz()>eR2CenterMax) s->GetElements()->Clear();
  }
}

DrawCorrMap()

void EdbViewMatch::DrawCorrMap

(

)

{
  bool batch = gROOT->IsBatch();
  if(eOutputFile) {
    Log(2,"DrawCorrMap","Save to file %s", eOutputFile->GetName());
    gROOT->SetBatch();
  }
  
  TCanvas *cc = new TCanvas("cc","Corrections map",1200,900);
      
  gStyle->SetPalette(1);
  gStyle->SetOptStat(0);
  gPad->SetGridx(1);
  gPad->SetGridy(1);
 
  float margin=10;
  double minXborder = eCorrMap.Xmin() - margin;
  double maxXborder = eCorrMap.Xmax() + margin;
  double minYborder = eCorrMap.Ymin() - margin;
  double maxYborder = eCorrMap.Ymax() + margin;
 
  TH2F *hh = new TH2F("hh","Corrections map",100,minXborder,maxXborder,100,minYborder,maxYborder);
  hh->GetXaxis()->SetTitle("X (#mum)");
  hh->GetYaxis()->SetTitle("Y (#mum)");
  hh->Draw();
 
  TBox *plate = new TBox(eCorrMap.Xmin(),eCorrMap.Ymin(),eCorrMap.Xmax(),eCorrMap.Ymax());
  plate->SetFillColor(16);
  plate->SetFillStyle(3001);
  plate->Draw();
 
  Double_t meanx=0, meany=0, wtot=0;
  float scale = 15.;
  int nc=eCorrMap.Ncell();
  for( int i=0; i<nc; i++ ) {
    int w = eCorrMap.Bin(i);    if(!w) continue;
    TArrayD *dxy = (TArrayD*)( eCorrMap.GetObject(i, 0 ) );
    float x = eCorrMap.Xj(i);
    float y = eCorrMap.Yj(i);
    float dx = (*dxy)[0];
    float dy = (*dxy)[1];
 
    TArrow *arrow = new TArrow(x,y,x+scale*dx,y+scale*dy,0.01);
    arrow->SetLineWidth(1);
    arrow->Draw();
    
    meanx += dx;
    meany += dy;
    wtot  += w;
  }
  
  meanx /= wtot;
  meany /= wtot;
  printf("\nmeanx = %g    meany = %g  wtot = %f\n", meanx, meany, wtot);
 
  eGMap.DrawH2("eGMap","entries/bin")->Write();
  eCorrMap.DrawH2("eCorrMap","entries/bin")->Write();
  cc->Write("corr_map");
  gROOT->SetBatch(batch);
}

DumpGr()

TNtuple * EdbViewMatch::DumpGr

(

const char *

name

)

{
  TNtuple *nt = new TNtuple( name,"grains dump","ig:n:xg:yg:rmin:nc:ic:xc:yc:xcv:ycv");
  int ngr = eGr.GetEntries();
  for(int ig=0; ig<ngr; ig++) {
    EdbSegment *s = ((EdbSegment*)eGr.UncheckedAt(ig));
    int nc = s->GetNelements();                                  if(!nc) continue;
    for( int ic=0; ic<nc; ic++ ) {
      EdbClMatch *c = (EdbClMatch *)(s->GetElements()->UncheckedAt(ic));
      nt->Fill( ig, s->GetNelements(), s->GetX0(), s->GetY0(), s->GetDz(), nc, ic, c->eX, c->eY, c->eXv, c->eYv );
    }
  }
  return nt;
}

GenerateCorrectionMatrix()

void EdbViewMatch::GenerateCorrectionMatrix

(

const char *

addfile

)

{
  int nentries = eCorrMap.Ncell();
  Double_t *vx   = new Double_t[nentries];
  Double_t *vy   = new Double_t[nentries];
  Double_t *vdx  = new Double_t[nentries];
  Double_t *vdy  = new Double_t[nentries];
 
  for(int i=0; i<nentries; i++) {
    int w = eCorrMap.Bin(i);    if(!w) continue;
    TArrayD *arr = (TArrayD*)eCorrMap.GetObject(i,0);
    if(!arr) { Log(1,"EdbViewMatch::GenerateCorrectionMatrix","ERROR: missed bin");  break; }
    vx[i]   = eCorrMap.Xj(i);
    vy[i]   = eCorrMap.Yj(i);
 
    vdx[i]  = arr->At(0);
    vdy[i]  = arr->At(1);
  }
 
  TGraph2D *gdx  = new TGraph2D("graphDX","graphDX",nentries,vx,vy,vdx);
  TGraph2D *gdy  = new TGraph2D("graphDY","graphDY",nentries,vx,vy,vdy);
  
  TH2F hdx("hdx","CorrectionMatrix dX",eNXpix, -eXpix*eNXpix/2., eXpix*eNXpix/2., eNYpix, -eYpix*eNYpix/2., eYpix*eNYpix/2.  );
  TH2F hdy("hdy","CorrectionMatrix dY",eNXpix, -eXpix*eNXpix/2., eXpix*eNXpix/2., eNYpix, -eYpix*eNYpix/2., eYpix*eNYpix/2.  );
  
  FILE *fadd=0;
  if(addfile) fadd = fopen(addfile,"r");
  char str[256];
  
   if(fadd) fgets(str, 256, fadd);    //skip first line
   int iadd, jadd;
   float xadd, yadd, dxadd, dyadd;
  
  FILE *fmatr = fopen("correction_matrix.txt","w");
  fprintf(fmatr,"%d %d  %f %f\n", eNXpix, eNYpix, eXpix, eYpix);
  
  for(int i=0; i<eNXpix; i++) {
    for(int j=0; j<eNYpix; j++) {
      float x = (i - 0.5*(eNXpix-1))*eXpix;
      float y = (j - 0.5*(eNYpix-1))*eYpix;
      //float dx = gdx->Interpolate(x,y);
      //float dy = gdy->Interpolate(x,y);
      float dx = ((TArrayD *)eCorrMap.GetObject(x,y,0))->At(0);
      float dy = ((TArrayD *)eCorrMap.GetObject(x,y,0))->At(1);
      
      if(fadd) { 
        if(fgets(str, 256, fadd)==NULL) Log(1,"GenerateCorrectionMatrix","ERROR: addfile is not correct");
        if( sscanf(str, "%d %d %f %f %f %f", &iadd, &jadd, &xadd, &yadd, &dxadd, &dyadd) != 6 )  {
          Log(1,"GenerateCorrectionMatrix","ERROR: addfile is not correct");
          break;
        }
        dx += dxadd;
        dy += dyadd;
      }
      fprintf(fmatr,"%5d %5d %12.6f %12.6f  %11.6f %11.6f\n", i, j, x,y,dx,dy);
 
      hdx.Fill(x,y,dx);
      hdy.Fill(x,y,dy);
    }
  }
  fclose(fmatr);
  if(fadd) fclose(fadd);
  
  bool batch = gROOT->IsBatch();
  if(eOutputFile) {
    Log(2,"GenerateCorrectionMatrix","Save to file %s", eOutputFile->GetName());
    gROOT->SetBatch();
  }
 
  //TCanvas c("c","CorrectionMatrix", 700, 1200 );
  //c.Divide(1,2);
  //c.cd(1);  hdx.Draw("colz");
  //c.cd(2);  hdy.Draw("colz");
  //c.Write();
  hdx.Write("hdx");
  hdy.Write("hdy");
  gdx->Write("gdx");
  gdy->Write("gdy");
  gROOT->SetBatch(batch);
  
}

InitCorrMap()

void EdbViewMatch::InitCorrMap

(

)

{
  // matrix of the corrections in each entry is TArrayD with 2 entries:(dx,dy)
  float  mi[2] = { -eNXpix*Abs(eXpix)/2., -eNYpix*Abs(eYpix)/2. };
  float  ma[2] = {  eNXpix*Abs(eXpix)/2.,  eNYpix*Abs(eYpix)/2. };
  int    n[2] = { int((ma[0]-mi[0])/eCorrectionMatrixStepX), int((ma[1]-mi[1])/eCorrectionMatrixStepY) };
  
  eCorrectionMatrixStepX = (ma[0]-mi[0])/n[0];
  eCorrectionMatrixStepY = (ma[1]-mi[1])/n[1];
  n[0] = int((ma[0]-mi[0]+1.)/eCorrectionMatrixStepX);
  n[1] = int((ma[1]-mi[1]+1.)/eCorrectionMatrixStepY);
  
  eCorrMap.InitCell(20, n, mi, ma);
  eCorrMap.PrintStat();
  int nc=eCorrMap.Ncell();
  for( int i=0; i<nc; i++ ) {
    TArrayD *a = new TArrayD(2);
    eCorrMap.AddObject(i, (TObject*)a );
    eCorrMap.SetBin(i,0);
  }
 
}

InitGMap()

void EdbViewMatch::InitGMap

(

)

{
  float mi[2] = {-500,-400}, ma[2]={1500,1100}; 
  int    n[2] = { int((ma[0]-mi[0])/10), int((ma[1]-mi[1])/10) };
  eGMap.InitCell(250, n, mi, ma);
}

MakeDistortionMap()

void EdbViewMatch::MakeDistortionMap	(	const char *	fname,
		TEnv &	env,
		const char *	addfile = 0
	)

{
  EdbRun run(fname);
  SetPar(cenv);
  InitGMap();
  InitCorrMap();
  int      n = run.GetEntries();
 
  EdbView *v = run.GetView();
  run.GetEntry(0,1);
  //float X0 = v->GetXview();
  //float Y0 = v->GetYview();
  float X0 = 0;
  float Y0 = 0;
 
  for(int i=0; i<n; i++) {
    run.GetEntry(i,1,1);
    int ncl = v->Nclusters();
    printf("%6d ", ncl);
 
    for(int ic=0; ic<ncl; ic++) {
      EdbCluster *c = v->GetCluster(ic);
      if(c->GetArea()>eAreaMin&&c->GetArea()<eAreaMax) {
        AddCluster( c->eX, c->eY, c->eZ, v->GetXview()-X0, v->GetYview()-Y0, v->GetViewID(), c->GetFrame() );
      }
    }
  }
  printf("\n");
 
  CalculateGrRef();
  
  eOutputFile = new TFile("map.root","RECREATE");
  if(eDumpGr) { TNtuple *ntgr = DumpGr("gr"); ntgr->Write(); }
  CutGrRef();
  CalculateCorr();
  DrawCorrMap();
  GenerateCorrectionMatrix(addfile);
  if(eDumpGr) { TNtuple *ntgr = DumpGr("grcut"); ntgr->Write(); }
  eOutputFile->Close();
  Print();
}

Print()

void EdbViewMatch::Print

(

)

{
  printf("\n-----------------------------------------------------------------------------------------------------\n");
  int ncl = eCl.GetEntries();
  int ngr = eGr.GetEntries();
  int icgr=0,  iccl=0;
  for(int i=0; i<ngr; i++) {
    int nc = ((EdbSegment*)eGr.At(i))->GetNelements();
    if(nc>=eNClMin) { icgr++;  iccl+=nc; }
    
    //if(nc>500)   printf("grain %d  with  %d  cl\n", i, nc);
  }
  printf("\n%d clusters processed\n",ncl);
  printf("\n%d grains found\n",ngr);
  printf("\n%d grains used (ncl > %d and closer then %f to center) \n",icgr, eNClMin, eR2CenterMax );
  printf("\n%d clusters used in good grains\n",iccl);
  printf("Matrix definition: %d %d  %f %f view size: %f %f\n", eNXpix, eNYpix, eXpix, eYpix, eNXpix*eXpix, eNYpix*eYpix );
 
  eGMap.PrintStat();
  eCorrMap.PrintStat();
  printf("-----------------------------------------------------------------------------------------------------\n");
}

SetPar()

void EdbViewMatch::SetPar

(

TEnv &

env

)

{
  eNClMin      = env.GetValue("viewdist.NClMin" , 20);
  eR2CenterMax = env.GetValue("viewdist.R2CenterMax" , 15.);
  eRmax        = env.GetValue("viewdist.Rmax" ,  1.);
 
  eXpix  = env.GetValue("viewdist.Xpix" , 0.30625);
  eYpix  = env.GetValue("viewdist.Ypix" , 0.30714);
  eNXpix = env.GetValue("viewdist.NXpix", 1280);
  eNYpix = env.GetValue("viewdist.NYpix", 1024);
  
  eDumpGr = env.GetValue("viewdist.DumpGr", 0);
  
  sscanf( env.GetValue("viewdist.ClusterAreaLimits", "0 90000000"), "%f %f", &eAreaMin, &eAreaMax);
  sscanf( env.GetValue("viewdist.ClusterVolumeLimits", "0 90000000"), "%f %f", &eVolumeMin, &eVolumeMax);
  
  printf("\n----------------------- Processing Parameters ---------------------------\n");
  printf("eNClMin\t %d\n",eNClMin);
  printf("eR2CenterMax\t %6.2f\n",eR2CenterMax);
  printf("eRmax\t %f6.2\n",eRmax);
  printf("Pixel:  %9.7f x %9.7f \n", eXpix, eYpix );
  printf("Matrix: %d x %d pixels, \t  %15.7f x %15.7f microns", eNXpix, eNYpix, eXpix*eNXpix, eYpix*eNYpix);
  printf("Clusters Area limits: %7.0f  %7.0f \n",  eAreaMin, eAreaMax );
  printf("Clusters Volume limits: %7.0f  %7.0f \n",  eVolumeMin, eVolumeMax );
  printf("-------------------------------------------------------------------------\n\n");
  
}

SetPixelSize()

void EdbViewMatch::SetPixelSize ( float  xpix, float  ypix  )

inline

{ eXpix=xpix; eYpix=ypix; }

Member Data Documentation

eAreaMax

float EdbViewMatch::eAreaMax

eAreaMin

float EdbViewMatch::eAreaMin

eCl

TClonesArray EdbViewMatch::eCl

private

eCorrectionMatrixStepX

float EdbViewMatch::eCorrectionMatrixStepX

private

eCorrectionMatrixStepY

float EdbViewMatch::eCorrectionMatrixStepY

private

eCorrMap

EdbCell2 EdbViewMatch::eCorrMap

private

eDumpGr

bool EdbViewMatch::eDumpGr

eGMap

EdbCell2 EdbViewMatch::eGMap

private

eGr

TClonesArray EdbViewMatch::eGr

private

eNClMin

int EdbViewMatch::eNClMin

private

eNXpix

int EdbViewMatch::eNXpix

private

eNYpix

int EdbViewMatch::eNYpix

private

eOutputFile

TFile* EdbViewMatch::eOutputFile

eR2CenterMax

float EdbViewMatch::eR2CenterMax

private

eRmax

float EdbViewMatch::eRmax

private

eVolumeMax

float EdbViewMatch::eVolumeMax

eVolumeMin

float EdbViewMatch::eVolumeMin

eXpix

float EdbViewMatch::eXpix

private

eYpix

float EdbViewMatch::eYpix

private

---

The documentation for this class was generated from the following files:

• /home/antonio/fedra_doxygen/src/libEMC/EdbViewMatch.h
• /home/antonio/fedra_doxygen/src/libEMC/EdbViewMatch.cxx