ConstructPLUT.c

Go to the documentation of this file.

/*
 *  Copyright (C) 2005 Badri Krishnan, Alicia Sintes, Bernd Machenschalk
 *
 *  This program is free software; you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation; either version 2 of the License, or
 *  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with with program; see the file COPYING. If not, write to the 
 *  Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, 
 *  MA  02111-1307  USA
 */

/** \file  ConstructPLUT.c
 * \ingroup pulsarHough
 * \author Sintes, A and Krishnan, B
 * \brief Core routines for constructing the Hough Look-Up-Tables
 * \date $Date: 2007/08/21 14:27:41 $
 * Revision: $Id: ConstructPLUT.c,v 1.28 2007/08/21 14:27:41 badri Exp $
 *
 * History:   Created by Sintes June 7, 2001
 *            Modified by Badri Krishnan Feb 2003
 *-----------------------------------------------------------------------
 *
\par Description 

This module is the core of the Hough transform. The LAL function 
LALHOUGHConstructPLUT()
constructs the look up tables that will be used to build the so-called
partial-Hough maps. Each look up table is valid for a given sky-patch, time, and
frequency range around a certain  \verb@f0@ value. The look up table contains
all the necessary information regarding the borders of the annuli clipped on
the \lq projected' two dimensional sky-patch.

The inputs are:  HOUGHPatchGrid   containing the grid patch
information. This is independent of the sky location of the
patch, And  HOUGHParamPLUT  with all the other parameters needed.

The output is the look up table  HOUGHptfLUT  


\par Uses
\code
PLUTInitialize(HOUGHptfLUT  *)
FillPLUT(HOUGHParamPLUT *, HOUGHptfLUT *, HOUGHPatchGrid *)
CheckLeftCircle(REAL8, REAL8, REAL8, INT4 *, INT4 *, INT4 *, 
                HOUGHPatchGrid *)
CheckRightCircle(REAL8, REAL8, REAL8, INT4 *, INT4 *, INT4 *, 
                 HOUGHPatchGrid *)
DrawRightCircle(REAL8, REAL8, REAL8, INT4, INT4, COORType *, 
                HOUGHPatchGrid *)
DrawLeftCircle(REAL8, REAL8, REAL8, INT4, INT4, COORType *, 
               HOUGHPatchGrid *)
CheckLineCase(REAL8, REAL8, REAL8, REAL8 *, INT4 *)
FindExactLine(REAL8, REAL8, REAL8 *, REAL8 *)
FindLine(REAL8, REAL8, REAL8, REAL8 *, REAL8 *)
CheckLineIntersection(REAL8, REAL8, REAL8, INT4 *, INT4 *, INT4 *, 
                      HOUGHPatchGrid *)
DrawLine(REAL8, REAL8, REAL8,INT4, INT4, COORType *, HOUGHPatchGrid *)
Fill1Column(INT4, INT4*, HOUGHptfLUT *, HOUGHPatchGrid *)
FillCaseN1(INT4, INT4, INT4, HOUGHptfLUT *, HOUGHPatchGrid *)
FillCaseN2(INT4, INT4, HOUGHptfLUT *, HOUGHPatchGrid *)
FillCaseN3(INT4, INT4, INT4, INT4 *, HOUGHptfLUT *, HOUGHPatchGrid *)
FillCaseN4(INT4, INT4, HOUGHptfLUT *, HOUGHPatchGrid *)
FillCaseN5(INT4, INT4, INT4, HOUGHptfLUT *)
FillCaseN6(INT4, INT4, INT4, HOUGHptfLUT *, HOUGHPatchGrid *)
FillCaseN7(INT4, INT4, HOUGHptfLUT *, HOUGHPatchGrid *)
FillCaseN8(INT4, INT4, HOUGHptfLUT *, HOUGHPatchGrid *)
Fill1ColumnAnor(INT4, HOUGHptfLUT *, HOUGHPatchGrid *)
FillCaseA1(INT4, INT4, INT4, HOUGHptfLUT *)
FillCaseA2(INT4, INT4, INT4, HOUGHptfLUT *)
FillCaseA3(INT4, INT4, INT4, HOUGHptfLUT *, HOUGHPatchGrid *)
InitialCircleCase(INT4 *,REAL8, REAL8, REAL8, REAL8 *, INT4 *, INT4 *, 
                  HOUGHptfLUT *, HOUGHPatchGrid *)
SecondCircleCase(INT4, INT4*, REAL8, REAL8, REAL8, INT4, REAL8 *,
                 INT4*,INT4 *,INT4*, HOUGHptfLUT *, HOUGHPatchGrid *)
FollowCircleCase(INT4,INT4 *,REAL8,REAL8,REAL8,REAL8,REAL8,INT4 *,
                 INT4 *,INT4 *, HOUGHptfLUT *, HOUGHPatchGrid *)
InitialLineCase(INT4 *, REAL8, REAL8, REAL8, INT4 *, HOUGHptfLUT *, 
                HOUGHPatchGrid *)
SecondLineCase(INT4, INT4 *, REAL8, REAL8, REAL8, INT4 *, HOUGHptfLUT *, 
               HOUGHPatchGrid *)
FollowLineCase(INT4, INT4 *, REAL8, REAL8, REAL8, REAL8, INT4, INT4 *, 
                           HOUGHptfLUT *, HOUGHPatchGrid *)
\endcode


*/


/*
 * 1.  An author and Id block
 */

/************************************ <lalVerbatim file="ConstructPLUTCV">
Author: Sintes, A. M., Krishnan, B.
$Id: ConstructPLUT.c,v 1.28 2007/08/21 14:27:41 badri Exp $
************************************* </lalVerbatim> */

/*
 * 2. Commented block with the documetation of this module
 */


/* <lalLaTeX>
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Module \texttt{ConstructPLUT.c}}
\label{ss:ConstructPLUT.c}
Construction of the look up table for generating partial Hough maps assuming the
use of the stereographic projection.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection*{Prototypes}
\vspace{0.1in}
\input{ConstructPLUTD}
\index{\verb&LALHOUGHConstructPLUT()&}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection*{Macros (used only internally)}

\begin{verbatim}
#define MAX(A, B)  (((A) < (B)) ? (B) : (A))
#define MIN(A, B)  (((A) < (B)) ? (A) : (B))
#define cot(A)  (1./tan(A))
#define rint(x) floor((x)+0.5)
\end{verbatim}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection*{Static function declarations}
\begin{verbatim}
static void PLUTInitialize(HOUGHptfLUT  *);
static void FillPLUT(HOUGHParamPLUT *, HOUGHptfLUT *, HOUGHPatchGrid *);
static void CheckLeftCircle(REAL8, REAL8, REAL8, INT4 *, INT4 *, INT4 *, 
                            HOUGHPatchGrid *);
static void CheckRightCircle(REAL8, REAL8, REAL8, INT4 *, INT4 *, INT4 *, 
                             HOUGHPatchGrid *);
static void DrawRightCircle(REAL8, REAL8, REAL8, INT4, INT4, COORType *, 
                            HOUGHPatchGrid *);
static void DrawLeftCircle(REAL8, REAL8, REAL8, INT4, INT4, COORType *, 
                           HOUGHPatchGrid *);
static void CheckLineCase(REAL8, REAL8, REAL8, REAL8 *, INT4 *);
static void FindExactLine(REAL8, REAL8, REAL8 *, REAL8 *);
static void FindLine(REAL8, REAL8, REAL8, REAL8 *, REAL8 *);
static void CheckLineIntersection(REAL8, REAL8, REAL8, INT4 *, INT4 *, INT4 *, 
                                  HOUGHPatchGrid *);
static void DrawLine(REAL8, REAL8, REAL8,INT4, INT4, COORType *, HOUGHPatchGrid *);
static void Fill1Column(INT4, INT4*, HOUGHptfLUT *, HOUGHPatchGrid *);
static void FillCaseN1(INT4, INT4, INT4, HOUGHptfLUT *, HOUGHPatchGrid *);
static void FillCaseN2(INT4, INT4, HOUGHptfLUT *, HOUGHPatchGrid *);
static void FillCaseN3(INT4, INT4, INT4, INT4 *, HOUGHptfLUT *, HOUGHPatchGrid *);
static void FillCaseN4(INT4, INT4, HOUGHptfLUT *, HOUGHPatchGrid *);
static void FillCaseN5(INT4, INT4, INT4, HOUGHptfLUT *);
static void FillCaseN6(INT4, INT4, INT4, HOUGHptfLUT *, HOUGHPatchGrid *);
static void FillCaseN7(INT4, INT4, HOUGHptfLUT *, HOUGHPatchGrid *);
static void FillCaseN8(INT4, INT4, HOUGHptfLUT *, HOUGHPatchGrid *);
static void Fill1ColumnAnor(INT4, HOUGHptfLUT *, HOUGHPatchGrid *);
static void FillCaseA1(INT4, INT4, INT4, HOUGHptfLUT *);
static void FillCaseA2(INT4, INT4, INT4, HOUGHptfLUT *);
static void FillCaseA3(INT4, INT4, INT4, HOUGHptfLUT *, HOUGHPatchGrid *);
static void InitialCircleCase(INT4 *,REAL8, REAL8, REAL8, REAL8 *, INT4 *, INT4 *, 
                              HOUGHptfLUT *, HOUGHPatchGrid *);
static void SecondCircleCase(INT4, INT4*, REAL8, REAL8, REAL8, INT4, REAL8 *,
                             INT4*,INT4 *,INT4*, HOUGHptfLUT *, HOUGHPatchGrid *);
static void FollowCircleCase(INT4,INT4 *,REAL8,REAL8,REAL8,REAL8,REAL8,INT4 *,
                             INT4 *,INT4 *, HOUGHptfLUT *, HOUGHPatchGrid *);
static void InitialLineCase(INT4 *, REAL8, REAL8, REAL8, INT4 *, HOUGHptfLUT *, 
                            HOUGHPatchGrid *);
static void SecondLineCase(INT4, INT4 *, REAL8, REAL8, REAL8, INT4 *, HOUGHptfLUT *, 
                           HOUGHPatchGrid *);
static void FollowLineCase(INT4, INT4 *, REAL8, REAL8, REAL8, REAL8, INT4, INT4 *, 
                           HOUGHptfLUT *, HOUGHPatchGrid *);
\end{verbatim}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection*{Description}

This module is the core of the Hough transform. The LAL function \verb&LALHOUGHConstructPLUT()&
constructs the look up tables that will be used to build the so-called
partial-Hough maps. Each look up table is valid for a given sky-patch, time, and
frequency range around a certain  \verb@f0@ value. The look up table contains
all the necessary information regarding the borders of the annuli clipped on
the \lq projected' two dimensional sky-patch.

The inputs are:  \verb@HOUGHPatchGrid  *patch@ containing the grid patch
information. This is independent of the sky location of the
patch. And  \verb@HOUGHParamPLUT  *par@ with all the other parameters needed.

The output is:  \verb@HOUGHptfLUT   *lut@. The fields are:

\begin{description}
\item[\texttt{lut-> timeIndex }]  Time index of the {\sc lut}. 
\item[\texttt{lut->f0Bin }]  Frequency bin for which it has been
constructed.
\item[\texttt{lut->deltaF }]  Frequency resolution: \texttt{df=1/TCOH}
\item[\texttt{lut->nFreqValid }]  \verb@ = PIXERR * f0Bin * VEPI / VTOT@:
Number of frequencies where the {\sc lut} is
valid.
\item[\texttt{lut->iniBin }]  First bin affecting the patch with respect to
\verb@f0@.
\item[\texttt{lut->nBin }]  Exact number of bins affecting the patch.
\item[\texttt{lut->maxNBins }] Maximum number of bins affecting the patch (for
memory allocation purposes), i.e. length of \texttt{lut->bin}.
\item[\texttt{lut->maxNorders }] Maximum number of borders affecting the patch (for
memory allocation purposes), i.e. length of \texttt{lut->border}.
\item[\texttt{lut->border} ]  The annulus borders.
\item[\texttt{lut->bin} ]  Bin to border correspondence.
\end{description}



More detailed documentation can be found in the source code itself.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection*{Uses}
%\begin{verbatim}
%LALZDestroyVector()
%\end{verbatim}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection*{Notes}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\vfill{\footnotesize\input{ConstructPLUTCV}}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
</lalLaTeX> */


/*
 * 3. Includes. These should appear in the following order: 
 *    a. Standard library includes
 *    b. LDAS includes
 *    c. LAL includes
 */


#include <lal/LUT.h>


/*
 * 4. Assignment of Id string using NRCSID()  
 */

NRCSID (CONSTRUCTPLUTC, "$Id: ConstructPLUT.c,v 1.28 2007/08/21 14:27:41 badri Exp $");



/*
 * 5.a) Constants, structures (used only internally in this module)
 */

/*
 * 5.b) Type declarations (used only internally)
 */

/*
 * 5.c) Macros (used only internally)
 */

#define MAX(A, B)  (((A) < (B)) ? (B) : (A))
#define MIN(A, B)  (((A) < (B)) ? (A) : (B))
#define cot(A)  (1./tan(A))
#define rint(x) floor((x)+0.5)


/*
 * 5.d) Extern global variable declarations (strongly discouraged)
 */


/*
 * 5.e) static Global variables
 */


/*
 * 5.f) Static function declarations
 */

static void PLUTInitialize(HOUGHptfLUT  *);
static void FillPLUT(HOUGHParamPLUT *, HOUGHptfLUT *, HOUGHPatchGrid *);

static void CheckLeftCircle(REAL8,REAL8,REAL8,INT4 *,INT4 *,INT4 *, 
                          HOUGHPatchGrid *);
static void CheckRightCircle(REAL8,REAL8,REAL8,INT4 *,INT4 *,INT4 *, 
                           HOUGHPatchGrid *);
static void DrawRightCircle(REAL8,REAL8,REAL8,INT4,INT4, COORType *, 
                          HOUGHPatchGrid *);
static void DrawLeftCircle(REAL8,REAL8,REAL8,INT4,INT4,COORType *, 
                         HOUGHPatchGrid *);
static void CheckLineCase(REAL8, REAL8, REAL8,REAL8 *, INT4 *);
/* static void FindExactLine(REAL8,REAL8,REAL8 *,REAL8 *); */
static void FindLine(REAL8,REAL8,REAL8,REAL8 *,REAL8 *);
static void CheckLineIntersection(REAL8,REAL8,REAL8,INT4 *,INT4 *,INT4 *, 
                                HOUGHPatchGrid *);
static void DrawLine(REAL8, REAL8, REAL8,INT4, INT4, COORType *, HOUGHPatchGrid *);
static void Fill1Column(INT4, INT4*, HOUGHptfLUT *, HOUGHPatchGrid *);
static void FillCaseN1(INT4, INT4, INT4, HOUGHptfLUT *, HOUGHPatchGrid *);
static void FillCaseN2(INT4, INT4, HOUGHptfLUT *, HOUGHPatchGrid *);
static void FillCaseN3(INT4, INT4, INT4, INT4 *, HOUGHptfLUT *, HOUGHPatchGrid *);
static void FillCaseN4(INT4, INT4, HOUGHptfLUT *, HOUGHPatchGrid *);
static void FillCaseN5(INT4, INT4, INT4, HOUGHptfLUT *);
static void FillCaseN6(INT4, INT4, INT4, HOUGHptfLUT *, HOUGHPatchGrid *);
static void FillCaseN7(INT4, INT4, HOUGHptfLUT *, HOUGHPatchGrid *);
static void FillCaseN8(INT4, INT4, HOUGHptfLUT *, HOUGHPatchGrid *);
static void Fill1ColumnAnor(INT4, HOUGHptfLUT *, HOUGHPatchGrid *);
static void FillCaseA1(INT4, INT4, INT4, HOUGHptfLUT *);
static void FillCaseA2(INT4, INT4, INT4, HOUGHptfLUT *);
static void FillCaseA3(INT4, INT4, INT4, HOUGHptfLUT *, HOUGHPatchGrid *);

static void InitialCircleCase(INT4 *,REAL8, REAL8, REAL8, REAL8 *, INT4 *, INT4 *, 
                            HOUGHptfLUT *, HOUGHPatchGrid *);
static void SecondCircleCase(INT4, INT4*,REAL8,REAL8, REAL8,INT4,REAL8*,
                           INT4*,INT4 *,INT4*, HOUGHptfLUT *, HOUGHPatchGrid *);
static void FollowCircleCase(INT4,INT4 *,REAL8,REAL8,REAL8,REAL8,REAL8,INT4 *,
                           INT4 *,INT4 *, HOUGHptfLUT *, HOUGHPatchGrid *);
static void InitialLineCase(INT4 *, REAL8, REAL8, REAL8, INT4 *, HOUGHptfLUT *, 
                          HOUGHPatchGrid *);
static void SecondLineCase(INT4, INT4 *, REAL8, REAL8, REAL8, INT4 *, HOUGHptfLUT *, 
                         HOUGHPatchGrid *);
static void FollowLineCase(INT4, INT4 *,REAL8, REAL8, REAL8,REAL8,INT4,INT4 *, 
                         HOUGHptfLUT *, HOUGHPatchGrid *);


/*
 * 5.g)  The functions that make up the guts of this module
 */


/* >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< */
/*  ***************************** <lalVerbatim file="ConstructPLUTD"> */
void LALHOUGHConstructPLUT(LALStatus       *status,
                           HOUGHptfLUT     *lut,
                           HOUGHPatchGrid  *patch,
                           HOUGHParamPLUT  *par)  
{ /*  ************************************************ </lalVerbatim> */
 
  INT8    f0Bin;
  
  /* --------------------------------------------- */
  INITSTATUS (status, " LALHOUGHConstructPLUT", CONSTRUCTPLUTC);
  /*  ATTATCHSTATUSPTR (status); */

  /*   Make sure the arguments are not NULL: */ 
  /* use aborts instead of asserts */
  if (lut == NULL) {
    ABORT( status, LUTH_ENULL, LUTH_MSGENULL); 
  }
  /* ASSERT (lut, status, LUTH_ENULL, LUTH_MSGENULL); */

  if (patch == NULL) {
    ABORT( status, LUTH_ENULL, LUTH_MSGENULL); 
  }
  /* ASSERT (patch, status, LUTH_ENULL, LUTH_MSGENULL); */

  if (par == NULL) {
    ABORT( status, LUTH_ENULL, LUTH_MSGENULL); 
  }
  /* ASSERT (par,  status, LUTH_ENULL, LUTH_MSGENULL); */
  
  if (  fabs((REAL4)par->deltaF - (REAL4)patch->deltaF) > 1.0e-6) {
    ABORT( status, LUTH_EVAL, LUTH_MSGEVAL); 
  }
  /* ASSERT (par->deltaF == patch->deltaF,  status, LUTH_EVAL, LUTH_MSGEVAL); */

  /* -------------------------------------------   */

  f0Bin =  par->f0Bin;

  lut->deltaF = par->deltaF;
  lut->f0Bin  = f0Bin;
  
  lut->nFreqValid = par->nFreqValid;
  /* lut->nFreqValid = PIXERR * f0Bin *VEPI/VTOT; */

  /* -------------------------------------------   */
  
  PLUTInitialize(lut);  
  FillPLUT(par, lut, patch);
 
  
  /* Make sure number of bins makes sense with the dimensions  */
  if (lut->nBin <= 0) {
    ABORT( status, LUTH_ESIZE, LUTH_MSGESIZE); 
  }
  /* ASSERT (lut->nBin >0, status, LUTH_ESIZE, LUTH_MSGESIZE); */

  if (lut->nBin > lut->maxNBins) {
    ABORT( status, LUTH_ESIZE, LUTH_MSGESIZE); 
  }
  /* ASSERT (lut->nBin <= lut->maxNBins, status, LUTH_ESIZE, LUTH_MSGESIZE); */

  
  /* -------------------------------------------   */
  /* DETATCHSTATUSPTR (status); */
  
  /* normal exit */
  RETURN (status);
}


/* >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< */

/* >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< */
/*        >>>>>>>>     FUNCTIONS  ONLY  USED INTERNALLY     <<<<<<    */
/* >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< */

/* >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< */




/***************************************************************/
/*     Subroutine to initialize the partial LUT                */
/*     Authors:  Sintes, A.M                                   */
/*               June 7, 2001                                  */
/***************************************************************/

static void  PLUTInitialize(HOUGHptfLUT  *lut){
  UINT4  i;
  UINT4 maxNBins, maxNBorders;
  
  maxNBins = lut->maxNBins;
  maxNBorders = lut->maxNBorders;
  
  for(i=0;i<maxNBins;++i){

    /* check lut->bin is not null */
    if ( lut->bin + i == NULL ) {
      fprintf(stderr," pointer lut->bin+%d is null [ConstructPLUT.c %d]\n", i, __LINE__);
    }

    lut->bin[i].leftB1  = 0;
    lut->bin[i].rightB1 = 0;
    lut->bin[i].leftB2  = 0;
    lut->bin[i].rightB2 = 0;
    lut->bin[i].piece1max = -1; 
    lut->bin[i].piece1min = 0;
    lut->bin[i].piece2max = -1;
    lut->bin[i].piece2min = 0;
  }
  
  for(i=0;i<maxNBorders;++i){

    /* check lut->border is not null */
    if ( lut->border + i == NULL ) {
      fprintf(stderr," pointer lut->border+%d is null [ConstructPLUT.c %d]\n", i, __LINE__);
    }

    lut->border[i].yUpper = -1;
    lut->border[i].yLower = 0;
  }
  return;
}

/* >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< */
/* >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< */
/* >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< */

/****************************************************************/
/*       subroutine to fill the  LUT                            */
/*                                                              */
/*     Authors:  Sintes, A.M                                    */
/*               March 16, 2001                                 */
/****************************************************************/

/* ******************* Some explanations *********************** */
/* The program should go like this: */

/* First call subroutine XXXXXXXXXX.c */
/*    to calculate xi(t) for a given f0 (and demodulation parameters). */
/*    Then rotate that vector to our new coordinates, if it was not  */
/*    constructed in that way. */
   
/*    We need here |xi|, alpha, delta[-pi/2,pi/2],  on the sphere. */
/*       or its cartesian coordiantes  */

/*       The master equation to solve is: */
/*           f(t)-f0 = xi*n - xi*N */
/*        or */
/*           cos( phi)= [f(t)-f0 +xi*N]/|xi| */

/*        where:       xi*N= -|xi|* sin(delta) */

/*        Given a bin in the peakgram with */
/*                  Deltaf= f(t)(center bin)- f0= k df  */
/*                         (df => frequency resolution =1/TCOH) */
/*        Then */
/*           [cos(phi)]max = min[ 1, k (df/|xi|) + (xi*N +df/2)/|xi| ] */
/*        and */
/*           [cos(phi)]min = max[-1, k (df/|xi|) + (xi*N -df/2)/|xi| ] */
/*       Note  */
/*         if  [cos(phi)]max = 1,  do  not increase k !!! */
/*         if  [cos(phi)]min =-1,  do  not decrease k !!! */
          
/*           [cos(phi)]max(k) = min[ 1, df/|xi|+[cos(phi)]max(k-1) ] */
/*           [cos(phi)]min(k) = max[ 1, df/|xi|+[cos(phi)]min(k-1) ]  */
/* ********************************************************************* */
/*     for a bin,peak or border: */
        
/*     if cos(phi)=+-1  => nothing to be drawn! */
/*     otherwise */
/*         obtain phi in the interval [0, LAL_PI] */
/*         if delta+phi= pi/2, or delta-phi=pi/2 (carefull, with +2*pi) */
/*                   draw line */
/*         otherwise draw a circle */
/* ******************************************************************* */
/* a pixel j , corresponds to x(j)= patch->deltaX*(0.5+j- patch->xSide/2.) */
/* the nearest center of a pixel is j = round[ x/patch->deltaX +patch->xSide/2.-0.5]  */
/* ******************************************************************* */
/* The way to convert k into the bin index is the following:    */
/*    for k>=0, binindex= k                                     */
/*    for k<0,  binindex= nBin+iniBin-1-k                       */
/*   f0 corresponds to k=0                                      */
/* This will be used when reading the peakgram!                 */ 
/* ******************************************************************* */

static void  FillPLUT(HOUGHParamPLUT  *par, HOUGHptfLUT  *lut,
                    HOUGHPatchGrid  *patch){

  /********************************************************/
  /*    variables that need to be calculated before       */
  /********************************************************/

  REAL8 cosDelta;   /* = df/|xi|                                             */
  REAL8 cosPhiMax0; /* max val of cosPhi of freq bin containing patch center */
  REAL8 cosPhiMin0; /* mix val of cosPhi of freq bin containing patch center */
  REAL8 alpha;      /* = xi.alpha in the rotated coordinates                 */
  REAL8 delta;      /* = xi.delta                                            */
  REAL8 epsilon;    /* = 8 * LINERR/(f0Bin * VTOT * patchSize^2)             */
                   

  /********************************************************/
  UINT4 maxNBins;
  UINT4 maxNBorders;
  
  INT4 lastBorder =0;  /* counter of the last build border */ 
  INT4 currentBin =0;  /* counter of the bin studied       */

  INT4 ifailPlus  = 1; /* =1 (ok, continue to next bin), =0 (stop) */
  INT4 ifailMinus = 1; /* =1 (ok, continue to previous bin), =0 (stop) */

  INT4 directionPlus=-1;  /* = +1,or -1 */
  INT4 directionPlusZero=-1;  /* = +1,or -1 */
  INT4 directionMinus; /* = +1,or -1 */

  INT4 pathology=1;   /* =1 (normal), =0 (anormal case) */
  INT4 lineCase;    /* =1 line Case, =0 circle case */
  INT4 nBinPos;
  
  /********************************************************/

  REAL8 rCritic;
  REAL8 lambda; 
  REAL8 rcOldPlus;
  REAL8 rcOldMinus;
  REAL8 cosPhiMax, cosPhiMin, phi;
  REAL8 ang1,ang2;
  REAL8 eps;

  /********************************************************/
  maxNBins = lut->maxNBins;
  maxNBorders = lut->maxNBorders;
  
  alpha = par->xi.alpha;
  delta = par->xi.delta;
  cosDelta = par->cosDelta;
  cosPhiMax0 = par->cosPhiMax0;
  cosPhiMin0 = par->cosPhiMin0; 
  epsilon = par->epsilon;

  /********************************************************/

  /* Copy value of offset */
  lut->offset = par->offset;
  
  /********************************************************/

  lambda = 2* delta -LAL_PI*0.5;
  rCritic = 2* cos(lambda) /(1 - sin(lambda) );

  /* Initializing variables to irrelevant values, 
     since these values should never be used */
  rcOldPlus = rCritic;
  rcOldMinus = rCritic;
  /********************************************************/
  /*  starting with the (central) bin corresponding to:   */
  /*            Delta_f(t) =0 (k=0), border cosPhiMax     */
  /********************************************************/ 

  /* cosPhiMax = MIN(1, cosPhiMax0); */
  cosPhiMax = cosPhiMax0;
  
  if(cosPhiMax > 0.99999999){ /* avoid points or small circles */
    ifailPlus = 0;            /* do not go to the next bin */
    directionPlus = -1;
  } else{
    
    phi  = acos(cosPhiMax);   /* in the interval (0,PI) */
    ang1 = delta + phi;
    ang2 = delta - phi;
    
    /* check for lines, or numerical lines! */
    CheckLineCase(epsilon, ang1, ang2, &eps, &lineCase);

    if( lineCase ){
        /* line case */
      InitialLineCase(&lastBorder,alpha,delta,eps, &ifailPlus,lut, patch);
      directionPlus = -1;
    } else{
      /* circle case */
      InitialCircleCase(&lastBorder,alpha, ang1, ang2,  
                        &rcOldPlus, &directionPlus, &ifailPlus,lut, patch);
    }    
  }
  
  directionPlusZero= directionPlus;

  /********************************************************/
  /* moving to the other bins                             */
  /********************************************************/

  /********************************************************/
  /* Plus direction: increasing cosPhiMax  */
  /********************************************************/



  while (ifailPlus){
   
    ++currentBin;
    pathology = 1;
    
    /* Some checks */
    /* #ifdef CHECKHOUGHINDEX */
    if (currentBin > maxNBins || lastBorder>= maxNBorders ){ 
      fprintf(stderr,"currentBin=%d not in range 1 to maxNBins=%d\n"
              "or lastborder=%d >= maxNBorders=%d [ConstructPLUT.c %d]\n",
              currentBin,maxNBins,lastBorder,maxNBorders, __LINE__);
    }
    /* #endif */

    cosPhiMax = cosPhiMax + cosDelta;
    /* or cosPhiMax = MIN(1,cosPhiMax + cosDelta ); */

    if( cosPhiMax > 0.99999999){ /* check appropiate value */
      ifailPlus = 0;
    } else {

      phi  = acos(cosPhiMax); /* it is in the interval (0,pi) */
      ang1 = delta + phi;
      ang2 = delta - phi;
      
      /* check for lines, or numerical lines! */
      CheckLineCase(epsilon, ang1, ang2, &eps, &lineCase);
      
      if( lineCase ){
        /* line case */
        FollowLineCase(currentBin, &lastBorder,alpha,delta,eps,
                       rcOldPlus, directionPlus, &ifailPlus,lut, patch);
      } else{
        /* circle case */
        FollowCircleCase(currentBin, &lastBorder,alpha, ang1, ang2,rCritic, 
                         rcOldPlus, &pathology, &directionPlus, 
                         &ifailPlus, lut, patch);
      }    
    }
    
    if(pathology){
      Fill1Column(currentBin, &lastBorder, lut, patch);
    }else{
      Fill1ColumnAnor(currentBin, lut, patch); 
    }
  }
  
  nBinPos = currentBin;
  
  
  /********************************************************/
  /*  starting with the (central) bin corresponding to:   */
  /*            Delta_f(t) =0 (k=0), border cosPhiMin     */
  /********************************************************/ 

  /* cosPhiMin = MAX(-1, cosPhiMin0); */
  cosPhiMin = cosPhiMin0;

  if(cosPhiMin < -0.99999999){ /* avoid points or small circles */
    ifailMinus = 0;            /* do not go to the next bin */
  } else{

    phi  = acos(cosPhiMin);   /* in the interval (0,PI) */
    ang1 = delta + phi;
    ang2 = delta - phi;
    
    /* check for lines, or numerical lines! */
    CheckLineCase(epsilon, ang1, ang2, &eps, &lineCase);
    
    if( lineCase ){
      /* line case */
      SecondLineCase(currentBin,&lastBorder,alpha,delta,eps, &ifailMinus,
                     lut, patch);
      directionMinus = -1;
      pathology = 0;  
    } else{
      /* circle case */
      pathology = 1; /* provisionally */
      SecondCircleCase(currentBin, &lastBorder,alpha, ang1, ang2, 
                       directionPlusZero, &rcOldMinus, 
                       &pathology, &directionMinus, 
                       &ifailMinus, lut, patch);
    }    
  }

  /********************************************************/
  /*  the way to identify initial pathologies:            */
  /*    initial or second case being a line !             */
  /*    or two circles, both with direction=-1            */
  /********************************************************/

  if(pathology){
    Fill1Column(0, &lastBorder,lut, patch);
  }else{
    Fill1ColumnAnor(0, lut, patch);
  }

  /********************************************************/
  /* moving to the other bins                             */
  /********************************************************/

  /********************************************************/
  /* Minus direction: decreasing cosPhiMin                */
  /********************************************************/

  while(ifailMinus){
   
    ++currentBin;
    pathology = 1;
    
    /* Some checks */
    /* #ifdef CHECKHOUGHINDEX */
    if (currentBin > maxNBins || lastBorder>= maxNBorders ){ 
      fprintf(stderr,"currentBin=%d not in range 1 to maxNBins=%d\n"
              "or lastborder=%d >= maxNBorders=%d [ConstructPLUT.c %d]\n",
              currentBin,maxNBins,lastBorder,maxNBorders, __LINE__);
    }
    /* #endif */
    

    cosPhiMin = cosPhiMin - cosDelta;

    if( cosPhiMin < -0.99999999){ /* check appropiate value */
      ifailMinus = 0;
    } else {

      phi  = acos(cosPhiMin); /* it is in the interval (0,pi) */
      ang1 = delta + phi;
      ang2 = delta - phi;
      
      /* check for lines, or numerical lines! */
      CheckLineCase(epsilon, ang1, ang2, &eps, &lineCase);
      
      if( lineCase ){
        /* line case */
        FollowLineCase(currentBin, &lastBorder,alpha,delta,eps,
                       rcOldMinus, directionMinus, &ifailMinus,lut, patch);
      } else{
        /* circle case */
        FollowCircleCase(currentBin, &lastBorder,alpha, ang1, ang2,rCritic, 
                         rcOldMinus, &pathology, &directionMinus, 
                         &ifailMinus, lut, patch);
      }    
    }
    
    if(pathology){
      Fill1Column(currentBin, &lastBorder, lut, patch);
    }else{
      Fill1ColumnAnor(currentBin, lut, patch); 
    }
  }
  
  /********************************************************/
  /* set iniBin,nBin  etc */
  /********************************************************/

  lut->nBin = currentBin + 1;
  lut->iniBin = nBinPos - currentBin;

  return;
}

/* end of the subroutine*/
/* >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< */
/* >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< */
/* >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< */


/****************************************************************/
/*                                                              */
/*     from initialLineCase.c          March 16, 2001           */
/*                                                              */
/*       program to find and identify bins and borders          */
/*       in the initial line case                               */
/*                                                              */
/*     Authors:  Sintes, A.M                                    */
/*                                                              */
/****************************************************************/


/****************************************************************/

static void InitialLineCase(INT4  *lastBorderP, REAL8 alpha, REAL8  delta, 
                     REAL8 eps, INT4 *ifailP, HOUGHptfLUT  *lut,
                     HOUGHPatchGrid  *patch){


  INT4 lastBorder;
  

  REAL8  xA,yA;
  INT4 yymin,yymax;
  REAL8 xRel, slope;
  INT4 noIn;   /* if no intersection occurs noIn=0 */

  /* some paranoid error checking */
  if (patch == NULL ) {
    fprintf(stderr, "pointer patch is null [ConstructPLUT.c %d]\n", __LINE__);
  }

  if (lut == NULL ) {
    fprintf(stderr, "pointer lut is null [ConstructPLUT.c %d]\n",__LINE__);
  }

  if (ifailP == NULL ) {
    fprintf(stderr, "pointer ifailP is null [ConstructPLUT.c %d]\n",__LINE__);
  }

  if (lastBorderP == NULL ) {
    fprintf(stderr, "pointer lastBorderP is null [ConstructPLUT.c %d]\n",__LINE__);
  }
  lastBorder = *lastBorderP;



 /************************************************/
  FindLine(alpha, delta, eps, &xA, &yA);
  CheckLineIntersection(alpha, xA, yA, &yymin, &yymax, &noIn, patch);

  if( noIn ==0 ){
    *ifailP = 0;    /* =1 (ok), =0 (stop) */
    return;
  }
  ++lastBorder;

  if (yymin < 0) {
    fprintf(stderr,"WARNING: Fixing yymin (%d -> 0) [ConstructPLUT.c %d]\n",
            yymin, __LINE__);
    yymin = 0;
  }
  if (yymax >= patch->ySide) {
    fprintf(stderr,"WARNING: Fixing yymax (%d -> %d) [ConstructPLUT.c %d]\n",
            yymax, patch->ySide-1, __LINE__);
    yymax = patch->ySide - 1;
  }

  lut->border[lastBorder].yUpper = yymax;
  lut->border[lastBorder].yLower = yymin;

  DrawLine(alpha, xA, yA, yymin, yymax,
           &lut->border[lastBorder].xPixel[0] , patch);

  /************************************************/

  if( (alpha == LAL_PI*0.5) || (alpha == LAL_PI*1.5) ||
      (alpha == -LAL_PI*0.5) ){       /* horizontal line */

    if( yA < 0 ){ /* convention */
      lut->bin[0].rightB1 = lastBorder;
      lut->bin[1].leftB1  = lastBorder;
      lut->border[lastBorder].yCenter = 0; /*or smaller*/
    } else {
      lut->bin[0].leftB2  = lastBorder;
      lut->bin[1].rightB2 = lastBorder;
      lut->border[lastBorder].yCenter = patch->ySide -1;/*or bigger */
    }
  } 
  else { /* non- horizontal line */
    xRel = xA + tan(alpha)*yA;

    if( (alpha == 0) || (alpha == LAL_PI) || (alpha == -LAL_PI) ){ 
      /* vertical line */
      slope = +1; /* arbitrary number, just to avoid overflow */
    } else {
      slope = - cot(alpha);
    } 

    if( xRel < 0 ){
      lut->bin[0].leftB2 = lastBorder;
      lut->bin[1].rightB2= lastBorder;

      if ( slope < 0 ){
        lut->border[lastBorder].yCenter = 0; /*or smaller*/
      } else {
        lut->border[lastBorder].yCenter = patch->ySide -1;/*or bigger */
      }

    } else {
      lut->bin[0].rightB1 = lastBorder;
      lut->bin[1].leftB1  = lastBorder;
      
      if ( slope > 0 ){
        lut->border[lastBorder].yCenter = 0; /*or smaller*/
      } else {
        lut->border[lastBorder].yCenter = patch->ySide -1;/*or bigger */
      }
    }
  } /* end non-horizontal line */


  /************************************************/

  *lastBorderP =  lastBorder;
 
 return;
}

/* >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< */
/* >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< */
/* >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< */

/****************************************************************/
/*                                                              */
/*     from secondLineCase.c          March 16, 2001            */
/*                                                              */
/*       program to find and identify bins and borders          */
/*       in the second line case                                */
/*                                                              */
/*     Authors:  Sintes, A.M                                    */
/*                                                              */
/****************************************************************/


/****************************************************************/

static void SecondLineCase(INT4 currentBin, INT4  *lastBorderP, 
                    REAL8 alpha, REAL8  delta, 
                     REAL8 eps, INT4 *ifailP, HOUGHptfLUT *lut,
                     HOUGHPatchGrid  *patch){

  /* we are for sure in a pathological case. Border names are
     changed accordingly */

  INT4 lastBorder;
  

  REAL8  xA,yA;
  INT4 yymin,yymax;
  REAL8 xRel, slope;
  INT4 noIn;   /* if no intersection occurs noIn=0 */


  /* some paranoid error checking */
  if (patch == NULL ) {
    fprintf(stderr, "pointer patch is null [ConstructPLUT.c %d]\n",__LINE__);
  }

  if (lut == NULL ) {
    fprintf(stderr, "pointer lut is null [ConstructPLUT.c %d]\n",__LINE__);
  }

  if</