Changeset 289 in Idarraga

Timestamp:

Apr 27, 2012, 12:02:03 AM (12 years ago)

Author:

idarraga

Message:

 

Location:

mafalda

Files:

• BlobsFinder/BlobsFinder.cpp (modified) (7 diffs)
• BlobsFinder/BlobsFinder.h (modified) (4 diffs)
• MPXAlgo/MediPixAlgo.cpp (modified) (1 diff)
• MPXAlgo/MediPixAlgo.h (modified) (2 diffs)

mafalda/BlobsFinder/BlobsFinder.cpp

@@ -12 +12 @@
 #include "MAFTools/MAFTools.h"
 
+#include "MPXAlgo/Highlighter.h"
 
 using namespace MSG;
@@ -50 +51 @@
         btype = type;
 }
-*/
+ */
 
 /** 
@@ -192 +193 @@
                                 if(StartBlob(colItr, rowItr))
                                 {
+                                        // Fill cluster description
+                                        FillClusterDescription();
+
                                         // calculate one blob's properties
                                         int retVal = CalculateProperties();
@@ -249 +253 @@
         oneBlob.blobContent.clear();
         oneBlob.blobContentSet.clear();
+        oneBlob.clusterDescription.clear();
 
         oneBlob.bP.nPixels = 0;
@@ -630 +635 @@
         oneBlob.bP.fitCut = fitCut;
 
+        // A rotation might be unreliable on a highly pixelized object.  Here's a different approach
+        //  to calculate the circle and elipse areas
+
+        // check if the object needs rotation
+        if(oneBlob.bP.chisquare_OverDof > m_chisquare_OverDof_rotation)
+        {
+
+                // Find the line passing through the extremes of the blob (using the linear regression is not quite right)
+                vector< pair<Float_t, Float_t> > limitPixels = MAFTools::FindLimitPixels(oneBlob);
+                double slope = ( limitPixels[1].second - limitPixels[0].second ) / ( limitPixels[1].first - limitPixels[0].first );
+                double cut = oneBlob.bP.geoCenter_y - slope*oneBlob.bP.geoCenter_x;
+
+                // Now calculate the perpendicular distance to this line from all points.
+                // I will take this as the b of the ellipse
+                // Then using the perpline I do the same to find a of the ellipse
+
+                double maxb = 0.;
+                double maxa = 0.;
+                double dist = 0., distcut = 0.;
+                double perpslope = -1./slope;
+                double perpcut = oneBlob.bP.geoCenter_y - perpslope*oneBlob.bP.geoCenter_x;
+                double cutdistance = MAFTools::CalcDistance(limitPixels[0].first, limitPixels[0].second, limitPixels[1].first, limitPixels[1].second) / oneBlob.bP.width_x;
+                double radiusSquared = 0.;
+
+                for (listItr = oneBlob.blobContent.begin() ; listItr != oneBlob.blobContent.end() ; listItr++) {
+
+                        // Tricky here ...
+                        // To determine b, I scan the points between the line at a perpendicular distance <= 1/5 of max pixels from the perp line
+                        // the opposite for a
+
+                        distcut = MAFTools::CalcPerpDistanceToLine(perpslope, perpcut, (*listItr).first );
+                        if(distcut <= cutdistance) {
+                                dist = MAFTools::CalcPerpDistanceToLine(slope, cut, (*listItr).first );
+                                if(dist > maxb) maxb = dist;
+                        }
+
+                        distcut = MAFTools::CalcPerpDistanceToLine(slope, cut, (*listItr).first );
+                        if(distcut <= cutdistance) {
+                                dist = MAFTools::CalcPerpDistanceToLine(perpslope, perpcut, (*listItr).first );
+                                if(dist > maxa) maxa = dist;
+                        }
+
+
+                        /*
+                        // Highlighters ... only for visualization
+                        // limits
+                        Highlighter * lineA = new Highlighter(limitPixels[0].first, limitPixels[0].second, "arrow", this);
+                        lineA->SetLineWidth(1);
+                        PullToStoreGateAccess(lineA, MPXDefs::DO_NOT_SERIALIZE_ME);
+
+                        Highlighter * lineB = new Highlighter(limitPixels[1].first, limitPixels[1].second, "arrow", this);
+                        lineB->SetLineWidth(1);
+                        PullToStoreGateAccess(lineB, MPXDefs::DO_NOT_SERIALIZE_ME);
+
+                        double x1 = oneBlob.bP.geoCenter_x - 5;
+                        double y1 = slope * x1 + cut;
+                        double x2 = oneBlob.bP.geoCenter_x + 5;
+                        double y2 = slope * x2 + cut;
+
+                        Highlighter * lineD = new Highlighter(x1, y1, x2, y2, "line", this);
+                        lineD->SetLineColor(kRed);
+                        lineD->SetLineWidth(1);
+                        PullToStoreGateAccess(lineD, MPXDefs::DO_NOT_SERIALIZE_ME);
+
+                        x1 = oneBlob.bP.geoCenter_x - 5;
+                        y1 = perpslope * x1 + perpcut;
+                        x2 = oneBlob.bP.geoCenter_x + 5;
+                        y2 = perpslope * x2 + perpcut;
+
+                        Highlighter * lineD2 = new Highlighter(x1, y1, x2, y2, "line", this);
+                        lineD2->SetLineColor(kGreen);
+                        lineD2->SetLineWidth(1);
+                        PullToStoreGateAccess(lineD2, MPXDefs::DO_NOT_SERIALIZE_ME);
+*/
+
+                }
+
+                oneBlob.bP.ellipseA = maxa;
+                oneBlob.bP.ellipseB = maxb;
+
+                // wx/2 is a, wy/2 is b.  AreaElipse = Pi*a*b
+                oneBlob.bP.ellipseArea = TMath::Pi()*(maxa*maxb)/4.0;
+
+                // fraction --> (sqrt(2)/2)
+                if(maxa >= maxb)
+                        radiusSquared = (2.*maxa*maxa) * (0.5/4.);
+                else
+                        radiusSquared = (2.*maxb*maxb) * (0.5/4.);
+
+                oneBlob.bP.circleArea = TMath::Pi()*radiusSquared;
+
+        }
+        else // no rotation needed
+        {
+                Double_t radiusSquared = (wx*wx + wy*wy) * 0.5/4.; // fraction --> (sqrt(2)/2)
+                oneBlob.bP.circleArea = TMath::Pi()*radiusSquared;
+                oneBlob.bP.ellipseArea = TMath::Pi()*(wx * wy)/4.0;
+                oneBlob.bP.ellipseA = wx;
+                oneBlob.bP.ellipseB = wy;
+        }
+
+        /*
         // check if the object needs rotation
         if(oneBlob.bP.chisquare_OverDof > m_chisquare_OverDof_rotation)
@@ -682 +789 @@
                 oneBlob.bP.ellipseB = wy;
         }
+         */
 
         // finally type initialized to _NOTYPE_ASSIGNED
@@ -1030 +1138 @@
 }
 
+void BlobsFinder::FillClusterDescription(){
+
+        list< pair < pair<Int_t, Int_t>, Int_t > > bc = oneBlob.GetBlobContent();
+        list< pair < pair<Int_t, Int_t>, Int_t > >::iterator i = bc.begin();
+
+        for ( ; i != bc.end() ; i++) {
+                oneBlob.clusterDescription.push_back( pair< pair<Int_t, Int_t>, Int_t > ( (*i).first, GetMatrixElement((*i).first) ) );
+        }
+
+}
+
 void BlobsFinder::CalcMinMaxGeoCenter(Float_t geox, Float_t geoy){
 

mafalda/BlobsFinder/BlobsFinder.h

@@ -144 +144 @@
   set< pair<Int_t, Int_t> > GetContentSet(){return blobContentSet;};
   list< pair < pair<Int_t, Int_t>, Int_t > > GetBlobContent(){return blobContent;};
+  list< pair < pair<Int_t, Int_t>, Int_t > > GetClusterDescription(){return clusterDescription;};
 
   // this will go private ... TODO !!!
@@ -151 +152 @@
   //                                               2  3  4
   // A <list> keeps entries in order
+  // WARNING The second member of the pair is not the TOT(or counts).  Is
+  //  the ordering in the cluster as defined above
   list< pair < pair<Int_t, Int_t>, Int_t > > blobContent;
   list< pair < pair<Int_t, Int_t>, Int_t > >::iterator blobItr;
@@ -156 +159 @@
   //  stdlib does better than me ;)
   set< pair<Int_t, Int_t> > blobContentSet;
+
+  // This list really contains the X,Y,C information
+  list< pair < pair<Int_t, Int_t>, Int_t > > clusterDescription;
 
   blobProperties bP;
@@ -221 +227 @@
   Float_t GetBalance(Int_t, Int_t, Int_t, Int_t, Float_t *, Float_t *, Int_t *, Int_t *, Float_t *, Float_t *);
   void CalcMinMaxGeoCenter(Float_t, Float_t);
+  void FillClusterDescription();
 
 private:

mafalda/MPXAlgo/MediPixAlgo.cpp

@@ -599 +599 @@
 }
 
+TCanvas * MediPixAlgo::DrawInSeparateWindow(vector<TGraph2D *> g){
+
+        ///////////////////////////////////////////////////////////////////
+        TCanvas * c11 = new TCanvas();
+        int ngraphs2D = (int)g.size();
+        int xdiv = floor( TMath::Sqrt(ngraphs2D) );
+        int ydiv = ceil( (double)ngraphs2D/(double)xdiv );
+        c11->Divide(xdiv, ydiv);
+
+        vector<TGraph2D *>::iterator itr = g.begin();
+        int cntr2 = 0;
+        for ( ; itr != g.end(); itr++) {
+                c11->cd(cntr2+1);
+                (*itr)->Draw("tri1");
+                cntr2++;
+        }
+
+        return c11;
+}
+
 #endif

mafalda/MPXAlgo/MediPixAlgo.h

@@ -18 +18 @@
 #include <TString.h>
 #include <TVector3.h>
+#include <TGraph2D.h>
+#include <TCanvas.h>
 
 #include "AnalysisCore/MediPixAnalysisCore.h"
@@ -193 +195 @@
   void ReadConfiguration();
 
+  // Extra drawing
+  TCanvas * DrawInSeparateWindow(vector<TGraph2D *>);
+
   /////////////////////////////
   // Report objects to the TBrowser