source: trunk/MagicSoft/Simulation/Detector/Camera/moments.cxx@ 6692

//=//////////////////////////////////////////////////////////////////////
//=
//= moments            
//=
//= @file        moments.cxx
//= @desc        Calculation of image parameters
//= @author      J C Gonzalez
//= @email       gonzalez@mppmu.mpg.de
//= @date        Thu May  7 16:24:22 1998
//=
//=----------------------------------------------------------------------
//=
//= Created: Thu May  7 16:24:22 1998
//= Author:  Jose Carlos Gonzalez
//= Purpose: Program for reflector simulation
//= Notes:   See files README for details
//=    
//=----------------------------------------------------------------------
//=
//= $RCSfile: moments.cxx,v $
//= $Revision: 1.1.1.1 $
//= $Author: harald $ 
//= $Date: 1999-11-05 11:59:33 $
//=
//=//////////////////////////////////////////////////////////////////////

// @T \newpage

//!@section Source code of |moments.cxx|.

/*!@{

  This section describes briefly the source code for the file
  |moments.cxx|. All the defines it uses are located in the file
  |moments.h|.

  @"*/

//!@subsection Includes and Global variables definition.

/*!@" 

  All the defines are located in the file |moments.h|.

  @"*/
 
//!@{

#include "moments.h"

//!@}

//!@subsection Definition of global variables.

//!@{

static int npix;                 //@< number of pixels
static float *q;                 //@< charges in the pixels
static float xm, ym;             //@< centroid (used in moments and lenwid)

//@: structure with information about the image
static Moments_Info m;           

//@: structure with information about islands
static Islands_Info is;          

//@: structure with information about lenwid
static LenWid_Info  lw;          

//!@}

//!@subsection The function |moments()|.

//!-----------------------------------------------------------
// @name moments               
//                                             
// @desc calculate moments on the image
//
// @var  n           Number of pixels
// @var  *image      Vector of ph.e.s in pixels
// @var  **pix       Array with information about the pixels
// @var  plateScale  Plate scale for the CT in use
// @var  flag        1: initialize; other: normal 
//
// @return           Pointer to structure Moments_Info
//
// @date Mon Sep 14 15:22:44 MET DST 1998
//------------------------------------------------------------
// @function

//!@{
Moments_Info * 
moments( int n, float *image, float **pix, 
         float plateScale, int flag )
{
  register int i, k;

  float x, y;
  float x2m, xym, y2m;
  float x3m, x2ym, xy2m, y3m;
  float zz, zd, zu, zv;
  float ax, ay, unitx, unity, sigmaax;
  float sx2, sxy, sy2;
  float sx3, sx2y, sxy2, sy3;
  
  if ( flag == 1 ) {
    q = new float[n];
    is.fi = new float[n];
    is.vislands = new float[n];
    is.islands = new int[n];
    is.isl = new int[n];
    for (i=1; i<n; ++i) {
      q[i] = is.fi[i] = is.vislands[i] = 0.0;
      is.islands[i] = is.isl[i] = 0;
    }
    return &m;
  } else {
    memcpy( q, image, sizeof(float) * n );
    /*
      for (i=1; i<n; ++i) 
      cout << q[i] << '\n';
      cout << endl << flush;
    */
  }

  // save number of pixels
  npix = n;
  
  // calculate sum of values
  xm = ym = 0.0;
  x2m = xym = y2m = 0.0;
  x3m = x2ym = xy2m = y3m = 0.0;
  m.charge = 0.0;
  
  for (i=0; i<npix; ++i)
    if ( q[i] > 0.0 ) {
      x = pix[i][0];
      y = pix[i][1];
      xm += x * q[i];
      ym += y * q[i];
      x2m += x * x * q[i];
      xym += x * y * q[i];
      y2m += y * y * q[i];
      x3m  += x * x * x * q[i];
      x2ym += x * x * y * q[i];
      xy2m += x * y * y * q[i];
      y3m  += y * y * y * q[i];
      m.charge += q[i];
    }

  xm *= plateScale;
  ym *= plateScale;
  x2m *= plateScale * plateScale;
  xym *= plateScale * plateScale;
  y2m *= plateScale * plateScale;
  x3m  *= plateScale * plateScale * plateScale;
  x2ym *= plateScale * plateScale * plateScale;
  xy2m *= plateScale * plateScale * plateScale;
  y3m  *= plateScale * plateScale * plateScale;
  
  //++++++++++++++++++++++++++++++++++++++++++++++++++ 
  // extremes and charges
  //--------------------------------------------------

  for (i=0; i<10; ++i) 
    m.maxs[i] = 0.0;

  for (i=0; i<npix; ++i) {
    if ( q[i] > m.maxs[0] ) {
      for (k=9; k>0; --k)
        m.maxs[k] = m.maxs[k-1];
      for (k=9; k>0; --k)
        m.nmaxs[k] = m.nmaxs[k-1];
      m.maxs[0] = q[i];
      m.nmaxs[0] = i;
    }
  }

  // calculates weighted position of the maximum (6 pixels)

  m.xmax = m.ymax = m.smax = 0.0;

  for (i=0; i<6; ++i) {
    m.xmax += pix[m.nmaxs[i]][0] * q[m.nmaxs[i]];
    m.ymax += pix[m.nmaxs[i]][1] * q[m.nmaxs[i]];
    m.smax += q[m.nmaxs[i]];
  }

  if (m.smax==0.) m.smax=1.;
  if (m.charge==0.) m.charge=1.;

  m.xmax = m.xmax * plateScale / m.smax;
  m.ymax = m.ymax * plateScale / m.smax;

  // calculate concentrations with 2,3,4...10 pixels

  m.conc[0] = q[ m.nmaxs[0] ] + q[ m.nmaxs[1] ];

  for (i=2; i<10; ++i) 
    m.conc[i-1] = m.conc[i-2] + q[ m.nmaxs[i] ];

  for (i=0; i<9; ++i) 
    m.conc[i] /= m.charge;  
  
  //++++++++++++++++++++++++++++++++++++++++++++++++++ 
  // 1st moments
  //--------------------------------------------------

  xm /= m.charge;
  ym /= m.charge;

  m.m1x = xm;
  m.m1y = ym;

  //++++++++++++++++++++++++++++++++++++++++++++++++++
  // 2nd moments
  //--------------------------------------------------

  x2m /= m.charge;
  xym /= m.charge;
  y2m /= m.charge;

  // around the origin
  m.m2xx = x2m;
  m.m2xy = xym;
  m.m2yy = y2m;
  
  // around the mean
  sx2 = x2m - SQR(xm);
  sxy = xym - xm * ym;
  sy2 = y2m - SQR(ym);

  m.m2cxx = sx2;
  m.m2cxy = sxy;
  m.m2cyy = sy2;

  //++++++++++++++++++++++++++++++++++++++++++++++++++
  // 3rd moments
  //--------------------------------------------------

  x3m  /= m.charge;
  x2ym /= m.charge;
  xy2m /= m.charge;
  y3m  /= m.charge;

  // around the origin
  m.m3xxx = x3m;
  m.m3xxy = x2ym;
  m.m3xyy = xy2m;
  m.m3yyy = y3m;

  // around the mean
  sx3  = x3m  - 3 * x2m * xm + 2 * xm * xm * xm;
  sx2y = x2ym - 2 * xym * xm + 2 * xm * xm * ym - x2m * ym;
  sxy2 = xy2m - 2 * xym * ym + 2 * xm * ym * ym - y2m * xm;
  sy3  = y3m  - 3 * y2m * ym + 2 * ym * ym * ym;

  m.m3cxxx = x3m;
  m.m3cxxy = x2ym;
  m.m3cxyy = xy2m;
  m.m3cyyy = y3m;

  //++++++++++++++++++++++++++++++++++++++++++++++++++ 
  // hillas' parameters
  //--------------------------------------------------

  zd = sy2 - sx2;
  zz = sqrt( SQR(zd) + 4.*SQR( sxy ));;
  if ( (zz < 1.e-6) || (sxy == 0.) ) {
    m.dist = -1.;
    return &m;
  }
  zu = 1.0 + zd / zz;
  zv = 2.0 - zu;

  /*
    a = (zd + zz) / (2 * sxy);
    b = ym - a * xm;
    
    m.length = sqrt( fabs( sx2 + 2 * a * sxy + a * a * sy2 ) / (1+a*a) );
    m.width  = sqrt( fabs( sx2 - 2 * a * sxy + a * a * sy2 ) / (1+a*a) );
    m.dist   = sqrt( SQR(xm) + SQR(ym) );
    m.xdist  = sqrt( SQR( m.xmax ) + SQR( m.ymax ) );
    m.azw    = sqrt( fabs( SQR(xm)* y2m - 2.* xm * ym * xym + x2m*SQR(ym) ) );
    m.miss   = fabs( b / sqrt(1+a*a) );
    m.alpha  = DEG( asin( m.miss / m.dist ) );
  */
  
  m.length = sqrt( fabs(sx2 + sy2 + zz) / 2. );
  m.width  = sqrt( fabs(sx2 + sy2 - zz) / 2. );
  m.dist   = sqrt( SQR(xm) + SQR(ym) );
  m.xdist  = sqrt( SQR(m.xmax) + SQR(m.ymax) );
  m.azw    = sqrt( fabs( SQR(xm)*y2m - 2.*xm*ym*xym + x2m*SQR(ym) ) ) / m.dist;
  m.miss   = sqrt( fabs( (SQR(xm)*zu + SQR(ym)*zv)/2. - (2.*sxy*xm*ym/zz) ) );
  m.alpha  = DEG( asin( m.miss/m.dist ) );
 

  /*
    length = sqrt( fabs(sx2 + sy2 + zz)  /2. );
    width  = sqrt( fabs(sx2 + sy2 - zz) / 2. );
    dist   = sqrt( SQR(xm) + SQR(ym) );
    xdist  = sqrt( SQR(m.xmax) + SQR(m.ymax) );
    azw    = sqrt( fabs( SQR(xm)*y2m - 2.*xm*ym*xym + x2m*SQR(ym) ) );
    miss   = sqrt( fabs( (SQR(xm)*zu + SQR(ym)*zv)/2. - (2.*sxy*xm*ym/zz) ) );
    alpha  = DEG( asin(miss/dist) );
  */  

  //++++++++++++++++++++++++++++++++++++++++++++++++++
  // asymetry
  //--------------------------------------------------
  
  unitx = sqrt(0.5*zv);
  unity = SGN( sxy )*sqrt(0.5*zu);

  if ( m.xdist > 0.0 ) {

    m.phi = acos((unitx*m.xmax + unity*m.ymax )/m.xdist);
  
    sigmaax = sx3*CUB(cos(m.phi)) + 3.0*sx2y*SQR(cos(m.phi))*sin(m.phi) +
      3.0*sxy2*cos(m.phi)*SQR(sin(m.phi)) + sy3*CUB(sin(m.phi));
    sigmaax = pow(fabs(sigmaax),0.3333333)*SGN(sigmaax);
    
    ax = sigmaax*unitx;
    ay = sigmaax*unity;
    m.asymx = (ax*m.xmax + ay*m.ymax)/(m.xdist*m.length*cos(m.phi));
    m.asymy = 0.0;

  } else {

    m.phi=-1000.0;
    m.asymx = -1000.0;
    m.asymy = -1000.0;

  }

  /*
        cout 
    << "length "<< length << endl 
    << "width  "<< width  << endl 
    << "dist   "<< dist   << endl 
    << "xdist  "<< xdist  << endl 
    << "azw    "<< azw    << endl 
    << "miss   "<< miss   << endl 
    << "alpha  "<< alpha  << endl << flush;
  */

  return &m;
}
//!@}

// @T \newpage

//!@subsection The function |islands()|.

//!-----------------------------------------------------------
// @name islands
//                                             
// @desc implementation of the "islands" algorithm
//
// @var  n             Number of pixels
// @var  f             Vector with the image (ph.e.s in pixels)
// @var  **pixneig     Array with indices of neighbour pixels
// @var  *npixneig     Vector with number of neighbours per pixel
// @var  cleanning     TRUE: remove spurious islands
// @var  ipixcut       Islands number cut
//
// @return           Pointer to structure Islands_Info
//
// @date Mon Sep 14 15:22:44 MET DST 1998
//------------------------------------------------------------
// @function

//!@{
Islands_Info *
islands( int n, float *f, int **pixneig, int *npixneig,
         int cleanning, int ipixcut)
{ 
  register int i;
  int j, k;
  int haschanged;
  
  is.numisl = 0;

  memcpy( is.fi, f, sizeof(float) * n );

  // be aware: here we use the side effect of ++
  // there are two possibilities of using the operator ++:
  // 1) a = ++i  => is.first increments i, then evaluates expresion
  // 2) a = i++  => is.first evaluates expresion, then increments i
  // we INTENTIONALLY use the second form

  // algorithm to isolate/detect is.islands
  j=1;
  for (i=0; i<n; ++i) 
    if ( is.fi[i]>0.0 ) {
      is.isl[i] = j;
      ++j;
    } else {
      is.isl[i] = 0;
    }
 
  haschanged = TRUE;
  while ( haschanged ) {
    haschanged = FALSE;
    for (i=0; i<n; ++i)
      if ( is.isl[i] > 0 )
        for (j=0; (j<npixneig[i]) && (pixneig[i][j]>-1); ++j)
          if ( (k=is.isl[pixneig[i][j]]) > 0 ) 
            if ( is.isl[i] > is.isl[k] ) {
              is.isl[i] = is.isl[k];
              haschanged=TRUE;
            } 
  }

  // count is.islands

  for (i=0;i<n;++i)
    is.islands[i] = 0;

  for (i=0;i<n;++i)
    is.vislands[i] = 0.0;

  for (i=0;i<n;++i)
    if (is.isl[i]>0) {
      is.islands[is.isl[i]]++;
      is.vislands[is.isl[i]] += is.fi[i];
    }
        
  for (i=0,j=0,is.numisl=0; i<n; ++i) {
    
    if (is.islands[i]>0) {
      j++;      
      //cout << '#' << j << ':' << is.islands[i] << "  q=" << is.vislands[i]
      //     << endl;
      
      if (is.islands[i] > ipixcut)
        is.numisl++;
    }
    
  }

  cout << j << '[' << is.numisl << "] is.islands\n" << flush;

  if ( cleanning ) {

    // cleanning image: pixcut = 3 (any is.island with <= 3 pixels is removed
    for (i=0;i<n;++i)
      if (is.islands[is.isl[i]] <= ipixcut)
        f[i] = 0.0;
  }

  return &is;
}
//!@}


//!@subsection The function |lenwid()|.

//!-----------------------------------------------------------
// @name lenwid
//                                             
// @desc calculation of extended length and width params.
//
// @var  n           Number of pixels
// @var  *image      Vector of ph.e.s in pixels
// @var  **pix       Array with information about the pixels
// @var  plateScale  Plate scale for the CT in use
// @var  flag        1: initialize; other: normal 
//
// @return           Pointer to structure LenWid_Info
//
// @date Mon Sep 14 15:22:44 MET DST 1998
//------------------------------------------------------------
// @function

//!@{
LenWid_Info * 
lenwid( int n, float *image, float **pix, 
        float plateScale,
        float max_distance)
{
  register int i, j, k;
  float chi, phi;
  float cp, sp;
  float px1[2], px2[2], py1[2], py2[2];
  float x1, x2, y1, y2;
  int sign_of_semiplane;
  float a, b, c;
  float dist_to_axis;
  float x, y;
  float wsum[4];
  float sum[4];
  float weight;
  float alpha;
  float radius, radius2;

  // calculate the radius of a circle with the same area of a pixel
  
  radius2 = max_distance*max_distance*cos(DEG30)*3.0 / M_PI;
  radius = sqrt(radius2);
  
  /* @comment
     We have now an image in the camera. In this image we have
     defined two axes, Xe and Ye. Given the definition of alpha,
     we define phi, which is the angle of the rotation that should
     be applied to the original axis X and Y to get, together with
     a translation to the point (xm, ym), the new axes Xe and Ye.
     @endcomment */

  chi = atan2(ym,xm);
  phi = m.alpha + chi;

  /* If the angle is phi, the rotation will be:
   *           / cos(phi)   sin(phi)\
   *  R(phi) = |                    |
   *           \-sin(phi)   cos(phi)/ 
   */

  cp=cos(phi);
  sp=sin(phi);

  /* The reference points for each axis will be px1,px2 and py1,py2
     We obtain these points by rotation and translation of the
     points [+-1000,0] and [0,+-1000] */

  /* Note! The rotation has to be R(-phi) */
  
  px1[0] = cp*1000 + xm;
  px1[1] = sp*1000 + ym;
    
  px2[0] = -cp*1000 + xm;
  px2[1] = -sp*1000 + ym;
    
  py1[0] = -sp*1000 + xm;
  py1[1] = cp*1000 + ym;
    
  py2[0] = sp*1000 + xm;
  py2[1] = -cp*1000 + ym;

  /* Now we have finally two points for each of the axes.
     We can now do, for each axis, and for each semi-plane
     it defines, the loop over the pixels */

  // Note that the possible values for sign_of_semiplane in the
  // next loops are precisely -1 and +1 

  for (i=0; i<4; ++i) {
    wsum[i] = sum[i] = 0.;
  }
  
  // first with the X, then with the Y

  for (k=1; k<=2; ++k) {
    
    if ( k == 1) {
      x1 = px1[0];
      y1 = px1[1];
      x2 = px2[0];
      y2 = px2[1];
    } else {
      x1 = py1[0];
      y1 = py1[1];
      x2 = py2[0];
      y2 = py2[1];
    }
  
    for ( sign_of_semiplane = -1;
          sign_of_semiplane < 2;
          sign_of_semiplane += 2 ) {
      
      // loop on pixels
      for ( i=0; i<n; ++i ) {
        
        // let's calculate the distance between the point and the axis
        
        x = pix[i][0];
        y = pix[i][1];
        
        a = (y2 - y1);
        b = (x1 - x2);
        c = (x1 * (y1-y2) + y1 * (x2 - x1));
      
        dist_to_axis = (a*x + b*y + c) / sqrt(a*a+b*b);
      
        // we have THREE cases:
      
        // (A)
        
        // if distance to the axis if larger than pixel diameter,
        // AND 
        // the semiplane is the WRONG one  -> forget that pixel
      
        if ( (fabs(dist_to_axis) > max_distance) &&
             (SGN(dist_to_axis) != sign_of_semiplane) )
          continue;
        
        // (B)
        
        // if distance to the axis if larger than pixel diameter,
        // AND 
        // the semiplane is the GOOD one  -> add this pixel
        
        if ( (fabs(dist_to_axis) > max_distance) &&
             (SGN(dist_to_axis) == sign_of_semiplane) ) {
          
          // here the sum
          
          weight = image[i];
          j = k+sign_of_semiplane;
          wsum[j] += weight * dist_to_axis * dist_to_axis;
          sum[j] += weight;
        
          continue;
        }
      
        // (C)
        // if we reach this point, that means that the center
        // of our pixel is too close to the axis, and we have
        // to feed it into the routine to check if the pixel
        // crosses the axis
      
        // ** NOTE ** NOTE ** NOTE ** NOTE ** NOTE ** NOTE ** NOTE
        // simplified algorithm
        // assume the pixels are circular, and takes
        // the fraction of the surface lying on the semiplane
        // ** NOTE ** NOTE ** NOTE ** NOTE ** NOTE ** NOTE ** NOTE
        
        // alpha
        alpha = 2*asin(sqrt(2*(radius-dist_to_axis)*radius -
                            radius2) / radius);
        
        // here the sum
        // the fraction is the fraction of the area inside the semiplane
        weight = image[i] * ( (alpha * radius2 / 2.0) / (M_PI * radius2));
        j = k+sign_of_semiplane;
        wsum[j] += weight * dist_to_axis * dist_to_axis;
        sum[j] += weight;   
      
      } // foreach pixel pixels
    
    } // foreach semiplane

  } // foreach axis

  lw.length1 = (sum[0] > 0.) ?  sqrt(wsum[0] / sum[0]) : -1;
  lw.width1  = (sum[1] > 0.) ?  sqrt(wsum[1] / sum[1]) : -1;
  lw.length2 = (sum[2] > 0.) ?  sqrt(wsum[2] / sum[2]) : -1;
  lw.width2  = (sum[3] > 0.) ?  sqrt(wsum[3] / sum[3]) : -1;
  
  lw.length1 *= plateScale;
  lw.width1  *= plateScale;
  lw.length2 *= plateScale;
  lw.width2  *= plateScale;

  return &lw;
}
//!@}


//!@subsection Auxiliary functions.

//!-----------------------------------------------------------
// @name crosspt
//                                             
// @desc calculate cross point of segments AB and CD
//
// @var ax         Coor. X of point A
// @var ay         Coor. Y of point A
// @var bx         Coor. X of point A
// @var by         Coor. Y of point A
// @var cx         Coor. X of point A
// @var cy         Coor. Y of point A
// @var dx         Coor. X of point A
// @var dy         Coor. Y of point A
// @var *pcrossx   Coor. X of cross point
// @var *pcrossy   Coor. Y of cross point
//
// @date Mon Mar  8 13:35:54 MET 1999
//------------------------------------------------------------
// @function

//!@{
void
crosspt( float ax, float ay,
         float bx, float by,
         float cx, float cy,
         float dx, float dy,
         float * pcrossx, float * pcrossy)
{
  float w, r;
  
  // the points A and B, and C and D define two segments (AB and CD)
  // the coordinates of these points are
  // A(ax,ay), B(bx,by), C(cx,cy), D(dx,dy)

  w=(bx-ax)*(dy-cy)-(by-ay)*(dx-cx);
  r=(ay-cy)*(dx-cx)-(ax-cx)*(dy-cy);

  *pcrossx = ax + r*(bx-ax)/w;
  *pcrossy = ay + r*(by-ay)/w;

}
//!@}

//=------------------------------------------------------------
//!@subsection Log of this file.

//!@{
//
// $Log: not supported by cvs2svn $
// Revision 1.2  1999/10/22 15:01:29  petry
// version sent to H.K. and N.M. on Fri Oct 22 1999
//
// Revision 1.1.1.1  1999/10/21 16:35:10  petry
// first synthesised version
//
// Revision 1.1  1999/03/08  10:04:06  gonzalez
// *** empty log message ***
//
// Revision 1.4  1999/03/02  09:56:14  gonzalez
// *** empty log message ***
//
// Revision 1.5  1999/03/15  14:59:10  gonzalez
// camera-1_1
//
//!@}

//=EOF