LALPulsar  6.1.0.1-89842e6
ConstructPLUT.c
Go to the documentation of this file.
1 /*
2  * Copyright (C) 2005 Badri Krishnan, Alicia Sintes, Bernd Machenschalk
3  *
4  * This program is free software; you can redistribute it and/or modify
5  * it under the terms of the GNU General Public License as published by
6  * the Free Software Foundation; either version 2 of the License, or
7  * (at your option) any later version.
8  *
9  * This program is distributed in the hope that it will be useful,
10  * but WITHOUT ANY WARRANTY; without even the implied warranty of
11  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12  * GNU General Public License for more details.
13  *
14  * You should have received a copy of the GNU General Public License
15  * along with with program; see the file COPYING. If not, write to the
16  * Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
17  * MA 02110-1301 USA
18  */
19 
20 /*
21  *
22  * History: Created by Sintes June 7, 2001
23  * Modified by Badri Krishnan Feb 2003
24  *-----------------------------------------------------------------------
25  */
26 
27 /**
28  * \file
29  * \ingroup LUT_h
30  * \author Sintes, A and Krishnan, B
31  * \brief Construction of the look up table for generating partial Hough maps assuming the
32  * use of the stereographic projection.
33  *
34  * ### Description ###
35  *
36  * This module is the core of the Hough transform. The LAL function LALHOUGHConstructPLUT()
37  * constructs the look up tables that will be used to build the so-called
38  * partial-Hough maps. Each look up table is valid for a given sky-patch, time, and
39  * frequency range around a certain \a f0 value. The look up table contains
40  * all the necessary information regarding the borders of the annuli clipped on
41  * the projected' two dimensional sky-patch.
42  *
43  * The inputs are: HOUGHPatchGrid containing the grid patch
44  * information. This is independent of the sky location of the
45  * patch, and HOUGHParamPLUT with all the other parameters needed.
46  *
47  * The output is the look up table HOUGHptfLUT
48  *
49  * More detailed documentation can be found in the source code itself.
50  *
51  */
52 
53 /*
54  * 3. Includes. These should appear in the following order:
55  * a. Standard library includes
56  * b. LDAS includes
57  * c. LAL includes
58  */
59 
60 
61 #include <lal/LUT.h>
62 
63 
64 /** \cond DONT_DOXYGEN */
65 
66 /*
67  * 5.a) Constants, structures (used only internally in this module)
68  */
69 
70 /*
71  * 5.b) Type declarations (used only internally)
72  */
73 
74 /*
75  * 5.c) Macros (used only internally)
76  */
77 
78 #define MAX(A, B) (((A) < (B)) ? (B) : (A))
79 #define MIN(A, B) (((A) < (B)) ? (A) : (B))
80 #define cot(A) (1./tan(A))
81 
82 /** \endcond */
83 
84 /*
85  * 5.d) Extern global variable declarations (strongly discouraged)
86  */
87 
88 
89 /*
90  * 5.e) static Global variables
91  */
92 
93 
94 /*
95  * 5.f) Static function declarations
96  */
97 
98 static void PLUTInitialize( HOUGHptfLUT * );
99 static void FillPLUT( HOUGHParamPLUT *, HOUGHptfLUT *, HOUGHPatchGrid * );
100 
101 static void CheckLeftCircle( REAL8, REAL8, REAL8, INT4 *, INT4 *, INT4 *,
102  HOUGHPatchGrid * );
103 static void CheckRightCircle( REAL8, REAL8, REAL8, INT4 *, INT4 *, INT4 *,
104  HOUGHPatchGrid * );
105 static void DrawRightCircle( REAL8, REAL8, REAL8, INT4, INT4, COORType *,
106  HOUGHPatchGrid * );
107 static void DrawLeftCircle( REAL8, REAL8, REAL8, INT4, INT4, COORType *,
108  HOUGHPatchGrid * );
109 static void CheckLineCase( REAL8, REAL8, REAL8, REAL8 *, INT4 * );
110 /* static void FindExactLine(REAL8,REAL8,REAL8 *,REAL8 *); */
111 static void FindLine( REAL8, REAL8, REAL8, REAL8 *, REAL8 * );
112 static void CheckLineIntersection( REAL8, REAL8, REAL8, INT4 *, INT4 *, INT4 *,
113  HOUGHPatchGrid * );
114 static void DrawLine( REAL8, REAL8, REAL8, INT4, INT4, COORType *, HOUGHPatchGrid * );
115 static void Fill1Column( INT4, UINT4 *, HOUGHptfLUT *, HOUGHPatchGrid * );
116 static void FillCaseN1( INT4, INT4, INT4, HOUGHptfLUT *, HOUGHPatchGrid * );
117 static void FillCaseN2( INT4, INT4, HOUGHptfLUT *, HOUGHPatchGrid * );
118 static void FillCaseN3( INT4, INT4, INT4, INT4 *, HOUGHptfLUT *, HOUGHPatchGrid * );
119 static void FillCaseN4( INT4, INT4, HOUGHptfLUT *, HOUGHPatchGrid * );
120 static void FillCaseN5( INT4, INT4, INT4, HOUGHptfLUT * );
121 static void FillCaseN6( INT4, INT4, INT4, HOUGHptfLUT *, HOUGHPatchGrid * );
122 static void FillCaseN7( INT4, INT4, HOUGHptfLUT *, HOUGHPatchGrid * );
123 static void FillCaseN8( INT4, INT4, HOUGHptfLUT *, HOUGHPatchGrid * );
124 static void Fill1ColumnAnor( INT4, HOUGHptfLUT *, HOUGHPatchGrid * );
125 static void FillCaseA1( INT4, INT4, INT4, HOUGHptfLUT * );
126 static void FillCaseA2( INT4, INT4, INT4, HOUGHptfLUT * );
127 static void FillCaseA3( INT4, INT4, INT4, HOUGHptfLUT *, HOUGHPatchGrid * );
128 
129 static void InitialCircleCase( UINT4 *, REAL8, REAL8, REAL8, REAL8 *, INT4 *, INT4 *,
130  HOUGHptfLUT *, HOUGHPatchGrid * );
131 static void SecondCircleCase( INT4, UINT4 *, REAL8, REAL8, REAL8, INT4, REAL8 *,
132  INT4 *, INT4 *, INT4 *, HOUGHptfLUT *, HOUGHPatchGrid * );
133 static void FollowCircleCase( INT4, UINT4 *, REAL8, REAL8, REAL8, REAL8, REAL8, INT4 *,
134  INT4 *, INT4 *, HOUGHptfLUT *, HOUGHPatchGrid * );
135 static void InitialLineCase( UINT4 *, REAL8, REAL8, REAL8, INT4 *, HOUGHptfLUT *,
136  HOUGHPatchGrid * );
137 static void SecondLineCase( INT4, UINT4 *, REAL8, REAL8, REAL8, INT4 *, HOUGHptfLUT *,
138  HOUGHPatchGrid * );
139 static void FollowLineCase( INT4, UINT4 *, REAL8, REAL8, REAL8, REAL8, INT4, INT4 *,
140  HOUGHptfLUT *, HOUGHPatchGrid * );
141 
142 
143 /*
144  * 5.g) The functions that make up the guts of this module
145  */
146 
147 
148 /* >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< */
149 void LALHOUGHConstructPLUT( LALStatus *status,
150  HOUGHptfLUT *lut,
151  HOUGHPatchGrid *patch,
152  HOUGHParamPLUT *par )
153 {
154 
155  INT8 f0Bin;
156 
157  /* --------------------------------------------- */
158  INITSTATUS( status );
159  /* ATTATCHSTATUSPTR (status); */
160 
161  /* Make sure the arguments are not NULL: */
162  /* use aborts instead of asserts */
163  if ( lut == NULL ) {
164  ABORT( status, LUTH_ENULL, LUTH_MSGENULL );
165  }
166  /* ASSERT (lut, status, LUTH_ENULL, LUTH_MSGENULL); */
167 
168  if ( lut->bin == NULL ) {
169  ABORT( status, LUTH_ENULL, LUTH_MSGENULL );
170  }
171 
172  if ( lut->border == NULL ) {
173  ABORT( status, LUTH_ENULL, LUTH_MSGENULL );
174  }
175 
176  if ( patch == NULL ) {
177  ABORT( status, LUTH_ENULL, LUTH_MSGENULL );
178  }
179  /* ASSERT (patch, status, LUTH_ENULL, LUTH_MSGENULL); */
180 
181  if ( par == NULL ) {
182  ABORT( status, LUTH_ENULL, LUTH_MSGENULL );
183  }
184  /* ASSERT (par, status, LUTH_ENULL, LUTH_MSGENULL); */
185 
186  if ( fabs( ( REAL4 )par->deltaF - ( REAL4 )patch->deltaF ) > 1.0e-6 ) {
187  ABORT( status, LUTH_EVAL, LUTH_MSGEVAL );
188  }
189  /* ASSERT (par->deltaF == patch->deltaF, status, LUTH_EVAL, LUTH_MSGEVAL); */
190 
191  /* ------------------------------------------- */
192 
193  f0Bin = par->f0Bin;
194 
195  lut->deltaF = par->deltaF;
196  lut->f0Bin = f0Bin;
197 
198  lut->nFreqValid = par->nFreqValid;
199  /* lut->nFreqValid = PIXERR * f0Bin *VEPI/VTOT; */
200 
201  /* ------------------------------------------- */
202 
203  PLUTInitialize( lut );
204  FillPLUT( par, lut, patch );
205 
206 
207  /* Make sure number of bins makes sense with the dimensions */
208  if ( lut->nBin <= 0 ) {
209  ABORT( status, LUTH_ESIZE, LUTH_MSGESIZE );
210  }
211  /* ASSERT (lut->nBin >0, status, LUTH_ESIZE, LUTH_MSGESIZE); */
212 
213  if ( lut->nBin > lut->maxNBins ) {
214  ABORT( status, LUTH_ESIZE, LUTH_MSGESIZE );
215  }
216  /* ASSERT (lut->nBin <= lut->maxNBins, status, LUTH_ESIZE, LUTH_MSGESIZE); */
217 
218 
219  /* ------------------------------------------- */
220  /* DETATCHSTATUSPTR (status); */
221 
222  /* normal exit */
223  RETURN( status );
224 }
225 
226 
227 /* >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< */
228 
229 /* >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< */
230 /* >>>>>>>> FUNCTIONS ONLY USED INTERNALLY <<<<<< */
231 /* >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< */
232 
233 /* >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< */
234 
235 
236 
237 
238 /***************************************************************/
239 /* Subroutine to initialize the partial LUT */
240 /* Authors: Sintes, A.M */
241 /* June 7, 2001 */
242 /***************************************************************/
243 
244 static void PLUTInitialize( HOUGHptfLUT *lut )
245 {
246  UINT4 i;
247  UINT4 maxNBins, maxNBorders;
248 
249  maxNBins = lut->maxNBins;
250  maxNBorders = lut->maxNBorders;
251 
252  for ( i = 0; i < maxNBins; ++i ) {
253  lut->bin[i].leftB1 = 0;
254  lut->bin[i].rightB1 = 0;
255  lut->bin[i].leftB2 = 0;
256  lut->bin[i].rightB2 = 0;
257  lut->bin[i].piece1max = -1;
258  lut->bin[i].piece1min = 0;
259  lut->bin[i].piece2max = -1;
260  lut->bin[i].piece2min = 0;
261  }
262 
263  for ( i = 0; i < maxNBorders; ++i ) {
264  lut->border[i].yUpper = -1;
265  lut->border[i].yLower = 0;
266  }
267  return;
268 }
269 
270 /* >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< */
271 /* >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< */
272 /* >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< */
273 
274 /****************************************************************/
275 /* subroutine to fill the LUT */
276 /* */
277 /* Authors: Sintes, A.M */
278 /* March 16, 2001 */
279 /****************************************************************/
280 
281 /* ******************* Some explanations *********************** */
282 /* The program should go like this: */
283 
284 /* First call subroutine XXXXXXXXXX.c */
285 /* to calculate xi(t) for a given f0 (and demodulation parameters). */
286 /* Then rotate that vector to our new coordinates, if it was not */
287 /* constructed in that way. */
288 
289 /* We need here |xi|, alpha, delta[-pi/2,pi/2], on the sphere. */
290 /* or its cartesian coordiantes */
291 
292 /* The master equation to solve is: */
293 /* f(t)-f0 = xi*n - xi*N */
294 /* or */
295 /* cos( phi)= [f(t)-f0 +xi*N]/|xi| */
296 
297 /* where: xi*N= -|xi|* sin(delta) */
298 
299 /* Given a bin in the peakgram with */
300 /* Deltaf= f(t)(center bin)- f0= k df */
301 /* (df => frequency resolution =1/TCOH) */
302 /* Then */
303 /* [cos(phi)]max = min[ 1, k (df/|xi|) + (xi*N +df/2)/|xi| ] */
304 /* and */
305 /* [cos(phi)]min = max[-1, k (df/|xi|) + (xi*N -df/2)/|xi| ] */
306 /* Note */
307 /* if [cos(phi)]max = 1, do not increase k !!! */
308 /* if [cos(phi)]min =-1, do not decrease k !!! */
309 
310 /* [cos(phi)]max(k) = min[ 1, df/|xi|+[cos(phi)]max(k-1) ] */
311 /* [cos(phi)]min(k) = max[ 1, df/|xi|+[cos(phi)]min(k-1) ] */
312 /* ********************************************************************* */
313 /* for a bin,peak or border: */
314 
315 /* if cos(phi)=+-1 => nothing to be drawn! */
316 /* otherwise */
317 /* obtain phi in the interval [0, LAL_PI] */
318 /* if delta+phi= pi/2, or delta-phi=pi/2 (carefull, with +2*pi) */
319 /* draw line */
320 /* otherwise draw a circle */
321 /* ******************************************************************* */
322 /* a pixel j , corresponds to x(j)= patch->deltaX*(0.5+j- patch->xSide/2.) */
323 /* the nearest center of a pixel is j = round[ x/patch->deltaX +patch->xSide/2.-0.5] */
324 /* ******************************************************************* */
325 /* The way to convert k into the bin index is the following: */
326 /* for k>=0, binindex= k */
327 /* for k<0, binindex= nBin+iniBin-1-k */
328 /* f0 corresponds to k=0 */
329 /* This will be used when reading the peakgram! */
330 /* ******************************************************************* */
331 
332 static void FillPLUT( HOUGHParamPLUT *par, HOUGHptfLUT *lut,
333  HOUGHPatchGrid *patch )
334 {
335 
336  /********************************************************/
337  /* variables that need to be calculated before */
338  /********************************************************/
339 
340  REAL8 cosDelta; /* = df/|xi| */
341  REAL8 cosPhiMax0; /* max val of cosPhi of freq bin containing patch center */
342  REAL8 cosPhiMin0; /* mix val of cosPhi of freq bin containing patch center */
343  REAL8 alpha; /* = xi.alpha in the rotated coordinates */
344  REAL8 delta; /* = xi.delta */
345  REAL8 epsilon; /* = 8 * LINERR/(f0Bin * VTOT * patchSize^2) */
346 
347 
348  /********************************************************/
349  UINT4 maxNBins;
350  UINT4 maxNBorders;
351 
352  UINT4 lastBorder = 0; /* counter of the last build border */
353  UINT4 currentBin = 0; /* counter of the bin studied */
354 
355  INT4 ifailPlus = 1; /* =1 (ok, continue to next bin), =0 (stop) */
356  INT4 ifailMinus = 1; /* =1 (ok, continue to previous bin), =0 (stop) */
357 
358  INT4 directionPlus = -1; /* = +1,or -1 */
359  INT4 directionPlusZero = -1; /* = +1,or -1 */
360  INT4 directionMinus; /* = +1,or -1 */
361 
362  INT4 pathology = 1; /* =1 (normal), =0 (anormal case) */
363  INT4 lineCase; /* =1 line Case, =0 circle case */
364  INT4 nBinPos;
365 
366  /********************************************************/
367 
368  REAL8 rCritic;
369  REAL8 lambda;
370  REAL8 rcOldPlus;
371  REAL8 rcOldMinus;
372  REAL8 cosPhiMax, cosPhiMin, phi;
373  REAL8 ang1, ang2;
374  REAL8 eps;
375 
376  /********************************************************/
377  maxNBins = lut->maxNBins;
378  maxNBorders = lut->maxNBorders;
379 
380  alpha = par->xi.alpha;
381  delta = par->xi.delta;
382  cosDelta = par->cosDelta;
383  cosPhiMax0 = par->cosPhiMax0;
384  cosPhiMin0 = par->cosPhiMin0;
385  epsilon = par->epsilon;
386 
387  /********************************************************/
388 
389  /* Copy value of offset */
390  lut->offset = par->offset;
391 
392  /********************************************************/
393 
394  lambda = 2 * delta - LAL_PI * 0.5;
395  rCritic = 2 * cos( lambda ) / ( 1 - sin( lambda ) );
396 
397  /* Initializing variables to irrelevant values,
398  since these values should never be used */
399  rcOldPlus = rCritic;
400  rcOldMinus = rCritic;
401  /********************************************************/
402  /* starting with the (central) bin corresponding to: */
403  /* Delta_f(t) =0 (k=0), border cosPhiMax */
404  /********************************************************/
405 
406  /* cosPhiMax = MIN(1, cosPhiMax0); */
407  cosPhiMax = cosPhiMax0;
408 
409  if ( cosPhiMax > 0.99999999 ) { /* avoid points or small circles */
410  ifailPlus = 0; /* do not go to the next bin */
411  directionPlus = -1;
412  } else {
413 
414  phi = acos( cosPhiMax ); /* in the interval (0,PI) */
415  ang1 = delta + phi;
416  ang2 = delta - phi;
417 
418  /* check for lines, or numerical lines! */
419  CheckLineCase( epsilon, ang1, ang2, &eps, &lineCase );
420 
421  if ( lineCase ) {
422  /* line case */
423  InitialLineCase( &lastBorder, alpha, delta, eps, &ifailPlus, lut, patch );
424  directionPlus = -1;
425  } else {
426  /* circle case */
427  InitialCircleCase( &lastBorder, alpha, ang1, ang2,
428  &rcOldPlus, &directionPlus, &ifailPlus, lut, patch );
429  }
430  }
431 
432  directionPlusZero = directionPlus;
433 
434  /********************************************************/
435  /* moving to the other bins */
436  /********************************************************/
437 
438  /********************************************************/
439  /* Plus direction: increasing cosPhiMax */
440  /********************************************************/
441 
442 
443 
444  while ( ifailPlus ) {
445 
446  ++currentBin;
447  pathology = 1;
448 
449  /* Some checks */
450  /* #ifdef CHECKHOUGHINDEX */
451  if ( currentBin > maxNBins || lastBorder >= maxNBorders ) {
452  fprintf( stderr, "currentBin=%d not in range 1 to maxNBins=%d\n"
453  "or lastborder=%d >= maxNBorders=%d [ConstructPLUT.c %d]\n",
454  currentBin, maxNBins, lastBorder, maxNBorders, __LINE__ );
455  }
456  /* #endif */
457 
458  cosPhiMax = cosPhiMax + cosDelta;
459  /* or cosPhiMax = MIN(1,cosPhiMax + cosDelta ); */
460 
461  if ( cosPhiMax > 0.99999999 ) { /* check appropiate value */
462  ifailPlus = 0;
463  } else {
464 
465  phi = acos( cosPhiMax ); /* it is in the interval (0,pi) */
466  ang1 = delta + phi;
467  ang2 = delta - phi;
468 
469  /* check for lines, or numerical lines! */
470  CheckLineCase( epsilon, ang1, ang2, &eps, &lineCase );
471 
472  if ( lineCase ) {
473  /* line case */
474  FollowLineCase( currentBin, &lastBorder, alpha, delta, eps,
475  rcOldPlus, directionPlus, &ifailPlus, lut, patch );
476  } else {
477  /* circle case */
478  FollowCircleCase( currentBin, &lastBorder, alpha, ang1, ang2, rCritic,
479  rcOldPlus, &pathology, &directionPlus,
480  &ifailPlus, lut, patch );
481  }
482  }
483 
484  if ( pathology ) {
485  Fill1Column( currentBin, &lastBorder, lut, patch );
486  } else {
487  Fill1ColumnAnor( currentBin, lut, patch );
488  }
489  }
490 
491  nBinPos = currentBin;
492 
493 
494  /********************************************************/
495  /* starting with the (central) bin corresponding to: */
496  /* Delta_f(t) =0 (k=0), border cosPhiMin */
497  /********************************************************/
498 
499  /* cosPhiMin = MAX(-1, cosPhiMin0); */
500  cosPhiMin = cosPhiMin0;
501 
502  if ( cosPhiMin < -0.99999999 ) { /* avoid points or small circles */
503  ifailMinus = 0; /* do not go to the next bin */
504  } else {
505 
506  phi = acos( cosPhiMin ); /* in the interval (0,PI) */
507  ang1 = delta + phi;
508  ang2 = delta - phi;
509 
510  /* check for lines, or numerical lines! */
511  CheckLineCase( epsilon, ang1, ang2, &eps, &lineCase );
512 
513  if ( lineCase ) {
514  /* line case */
515  SecondLineCase( currentBin, &lastBorder, alpha, delta, eps, &ifailMinus,
516  lut, patch );
517  directionMinus = -1;
518  pathology = 0;
519  } else {
520  /* circle case */
521  pathology = 1; /* provisionally */
522  SecondCircleCase( currentBin, &lastBorder, alpha, ang1, ang2,
523  directionPlusZero, &rcOldMinus,
524  &pathology, &directionMinus,
525  &ifailMinus, lut, patch );
526  }
527  }
528 
529  /********************************************************/
530  /* the way to identify initial pathologies: */
531  /* initial or second case being a line ! */
532  /* or two circles, both with direction=-1 */
533  /********************************************************/
534 
535  if ( pathology ) {
536  Fill1Column( 0, &lastBorder, lut, patch );
537  } else {
538  Fill1ColumnAnor( 0, lut, patch );
539  }
540 
541  /********************************************************/
542  /* moving to the other bins */
543  /********************************************************/
544 
545  /********************************************************/
546  /* Minus direction: decreasing cosPhiMin */
547  /********************************************************/
548 
549  while ( ifailMinus ) {
550 
551  ++currentBin;
552  pathology = 1;
553 
554  /* Some checks */
555  /* #ifdef CHECKHOUGHINDEX */
556  if ( currentBin > maxNBins || lastBorder >= maxNBorders ) {
557  fprintf( stderr, "currentBin=%d not in range 1 to maxNBins=%d\n"
558  "or lastborder=%d >= maxNBorders=%d [ConstructPLUT.c %d]\n",
559  currentBin, maxNBins, lastBorder, maxNBorders, __LINE__ );
560  }
561  /* #endif */
562 
563 
564  cosPhiMin = cosPhiMin - cosDelta;
565 
566  if ( cosPhiMin < -0.99999999 ) { /* check appropiate value */
567  ifailMinus = 0;
568  } else {
569 
570  phi = acos( cosPhiMin ); /* it is in the interval (0,pi) */
571  ang1 = delta + phi;
572  ang2 = delta - phi;
573 
574  /* check for lines, or numerical lines! */
575  CheckLineCase( epsilon, ang1, ang2, &eps, &lineCase );
576 
577  if ( lineCase ) {
578  /* line case */
579  FollowLineCase( currentBin, &lastBorder, alpha, delta, eps,
580  rcOldMinus, directionMinus, &ifailMinus, lut, patch );
581  } else {
582  /* circle case */
583  FollowCircleCase( currentBin, &lastBorder, alpha, ang1, ang2, rCritic,
584  rcOldMinus, &pathology, &directionMinus,
585  &ifailMinus, lut, patch );
586  }
587  }
588 
589  if ( pathology ) {
590  Fill1Column( currentBin, &lastBorder, lut, patch );
591  } else {
592  Fill1ColumnAnor( currentBin, lut, patch );
593  }
594  }
595 
596  /********************************************************/
597  /* set iniBin,nBin etc */
598  /********************************************************/
599 
600  lut->nBin = currentBin + 1;
601  lut->iniBin = nBinPos - currentBin;
602 
603  return;
604 }
605 
606 /* end of the subroutine*/
607 /* >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< */
608 /* >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< */
609 /* >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< */
610 
611 
612 /****************************************************************/
613 /* */
614 /* from initialLineCase.c March 16, 2001 */
615 /* */
616 /* program to find and identify bins and borders */
617 /* in the initial line case */
618 /* */
619 /* Authors: Sintes, A.M */
620 /* */
621 /****************************************************************/
622 
623 
624 /****************************************************************/
625 
626 static void InitialLineCase( UINT4 *lastBorderP, REAL8 alpha, REAL8 delta,
627  REAL8 eps, INT4 *ifailP, HOUGHptfLUT *lut,
628  HOUGHPatchGrid *patch )
629 {
630 
631 
632  INT4 lastBorder;
633 
634 
635  REAL8 xA, yA;
636  INT4 yymin, yymax;
637  REAL8 xRel, slope;
638  INT4 noIn; /* if no intersection occurs noIn=0 */
639 
640  /* some paranoid error checking */
641  if ( patch == NULL ) {
642  fprintf( stderr, "pointer patch is null [ConstructPLUT.c %d]\n", __LINE__ );
643  }
644 
645  if ( lut == NULL ) {
646  fprintf( stderr, "pointer lut is null [ConstructPLUT.c %d]\n", __LINE__ );
647  }
648 
649  if ( ifailP == NULL ) {
650  fprintf( stderr, "pointer ifailP is null [ConstructPLUT.c %d]\n", __LINE__ );
651  }
652 
653  if ( lastBorderP == NULL ) {
654  fprintf( stderr, "pointer lastBorderP is null [ConstructPLUT.c %d]\n", __LINE__ );
655  }
656  lastBorder = *lastBorderP;
657 
658 
659 
660  /************************************************/
661  FindLine( alpha, delta, eps, &xA, &yA );
662  CheckLineIntersection( alpha, xA, yA, &yymin, &yymax, &noIn, patch );
663 
664  if ( noIn == 0 ) {
665  *ifailP = 0; /* =1 (ok), =0 (stop) */
666  return;
667  }
668  ++lastBorder;
669 
670  if ( yymin < 0 ) {
671  fprintf( stderr, "WARNING: Fixing yymin (%d -> 0) [ConstructPLUT.c %d]\n",
672  yymin, __LINE__ );
673  yymin = 0;
674  }
675  if ( yymax >= patch->ySide ) {
676  fprintf( stderr, "WARNING: Fixing yymax (%d -> %d) [ConstructPLUT.c %d]\n",
677  yymax, patch->ySide - 1, __LINE__ );
678  yymax = patch->ySide - 1;
679  }
680 
681  lut->border[lastBorder].yUpper = yymax;
682  lut->border[lastBorder].yLower = yymin;
683 
684  DrawLine( alpha, xA, yA, yymin, yymax,
685  &lut->border[lastBorder].xPixel[0], patch );
686 
687  /************************************************/
688 
689  if ( ( alpha == LAL_PI * 0.5 ) || ( alpha == LAL_PI * 1.5 ) ||
690  ( alpha == -LAL_PI * 0.5 ) ) { /* horizontal line */
691 
692  if ( yA < 0 ) { /* convention */
693  lut->bin[0].rightB1 = lastBorder;
694  lut->bin[1].leftB1 = lastBorder;
695  lut->border[lastBorder].yCenter = 0; /*or smaller*/
696  } else {
697  lut->bin[0].leftB2 = lastBorder;
698  lut->bin[1].rightB2 = lastBorder;
699  lut->border[lastBorder].yCenter = patch->ySide - 1; /*or bigger */
700  }
701  } else { /* non- horizontal line */
702  xRel = xA + tan( alpha ) * yA;
703 
704  if ( ( alpha == 0 ) || ( alpha == LAL_PI ) || ( alpha == -LAL_PI ) ) {
705  /* vertical line */
706  slope = +1; /* arbitrary number, just to avoid overflow */
707  } else {
708  slope = - cot( alpha );
709  }
710 
711  if ( xRel < 0 ) {
712  lut->bin[0].leftB2 = lastBorder;
713  lut->bin[1].rightB2 = lastBorder;
714 
715  if ( slope < 0 ) {
716  lut->border[lastBorder].yCenter = 0; /*or smaller*/
717  } else {
718  lut->border[lastBorder].yCenter = patch->ySide - 1; /*or bigger */
719  }
720 
721  } else {
722  lut->bin[0].rightB1 = lastBorder;
723  lut->bin[1].leftB1 = lastBorder;
724 
725  if ( slope > 0 ) {
726  lut->border[lastBorder].yCenter = 0; /*or smaller*/
727  } else {
728  lut->border[lastBorder].yCenter = patch->ySide - 1; /*or bigger */
729  }
730  }
731  } /* end non-horizontal line */
732 
733 
734  /************************************************/
735 
736  *lastBorderP = lastBorder;
737 
738  return;
739 }
740 
741 /* >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< */
742 /* >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< */
743 /* >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< */
744 
745 /****************************************************************/
746 /* */
747 /* from secondLineCase.c March 16, 2001 */
748 /* */
749 /* program to find and identify bins and borders */
750 /* in the second line case */
751 /* */
752 /* Authors: Sintes, A.M */
753 /* */
754 /****************************************************************/
755 
756 
757 /****************************************************************/
758 
759 static void SecondLineCase( INT4 currentBin, UINT4 *lastBorderP,
760  REAL8 alpha, REAL8 delta,
761  REAL8 eps, INT4 *ifailP, HOUGHptfLUT *lut,
762  HOUGHPatchGrid *patch )
763 {
764 
765  /* we are for sure in a pathological case. Border names are
766  changed accordingly */
767 
768  INT4 lastBorder;
769 
770 
771  REAL8 xA, yA;
772  INT4 yymin, yymax;
773  REAL8 xRel, slope;
774  INT4 noIn; /* if no intersection occurs noIn=0 */
775 
776 
777  /* some paranoid error checking */
778  if ( patch == NULL ) {
779  fprintf( stderr, "pointer patch is null [ConstructPLUT.c %d]\n", __LINE__ );
780  }
781 
782  if ( lut == NULL ) {
783  fprintf( stderr, "pointer lut is null [ConstructPLUT.c %d]\n", __LINE__ );
784  }
785 
786  if ( ifailP == NULL ) {
787  fprintf( stderr, "pointer ifailP is null [ConstructPLUT.c %d]\n", __LINE__ );
788  }
789 
790  if ( lastBorderP == NULL ) {
791  fprintf( stderr, "pointer lastBorderP is null [ConstructPLUT.c %d]\n", __LINE__ );
792  }
793  lastBorder = *lastBorderP;
794 
795  /************************************************/
796  FindLine( alpha, delta, eps, &xA, &yA );
797  CheckLineIntersection( alpha, xA, yA, &yymin, &yymax, &noIn, patch );
798 
799  if ( noIn == 0 ) {
800  *ifailP = 0; /* =1 (ok), =0 (stop) */
801  return;
802  }
803  ++lastBorder;
804 
805  if ( yymin < 0 ) {
806  fprintf( stderr, "WARNING: Fixing yymin (%d -> 0) [ConstructPLUT.c %d]\n",
807  yymin, __LINE__ );
808  yymin = 0;
809  }
810  if ( yymax >= patch->ySide ) {
811  fprintf( stderr, "WARNING: Fixing yymax (%d -> %d) [ConstructPLUT.c %d]\n",
812  yymax, patch->ySide - 1, __LINE__ );
813  yymax = patch->ySide - 1;
814  }
815 
816  lut->border[lastBorder].yUpper = yymax;
817  lut->border[lastBorder].yLower = yymin;
818 
819  DrawLine( alpha, xA, yA, yymin, yymax,
820  &lut->border[lastBorder].xPixel[0], patch );
821 
822  /************************************************/
823  /* all are pathological cases. The code is similar to the
824  InitialLineCase, except for the modifications:
825  rightB1 -> rightB2, leftB2 -> leftB1, where marked */
826 
827  if ( ( alpha == LAL_PI * 0.5 ) || ( alpha == LAL_PI * 1.5 ) ||
828  ( alpha == -LAL_PI * 0.5 ) ) { /* horizontal line */
829 
830  if ( yA < 0 ) { /* convention */
831  lut->bin[0].rightB2 = lastBorder; /* modified */
832  lut->bin[currentBin + 1].leftB1 = lastBorder;
833  lut->border[lastBorder].yCenter = 0; /*or smaller*/
834  } else {
835  lut->bin[0].leftB1 = lastBorder; /* modified */
836  lut->bin[currentBin + 1].rightB2 = lastBorder;
837  lut->border[lastBorder].yCenter = patch->ySide - 1; /*or bigger */
838  }
839  } else { /* non- horizontal line */
840  xRel = xA + tan( alpha ) * yA;
841 
842  if ( ( alpha == 0 ) || ( alpha == LAL_PI ) || ( alpha == -LAL_PI ) ) {
843  /* vertical line */
844  slope = +1; /* arbitrary number, just to avoid overflow */
845  } else {
846  slope = - cot( alpha );
847  }
848 
849  if ( xRel < 0 ) {
850  lut->bin[0].leftB1 = lastBorder; /* modified */
851  lut->bin[currentBin + 1].rightB2 = lastBorder;
852 
853  if ( slope < 0 ) {
854  lut->border[lastBorder].yCenter = 0; /*or smaller*/
855  } else {
856  lut->border[lastBorder].yCenter = patch->ySide - 1; /*or bigger */
857  }
858 
859  } else {
860  lut->bin[0].rightB2 = lastBorder; /* modified */
861  lut->bin[currentBin + 1].leftB1 = lastBorder;
862 
863  if ( slope > 0 ) {
864  lut->border[lastBorder].yCenter = 0; /*or smaller*/
865  } else {
866  lut->border[lastBorder].yCenter = patch->ySide - 1; /*or bigger */
867  }
868  }
869  } /* end non-horizontal line */
870 
871 
872  /************************************************/
873 
874  *lastBorderP = lastBorder;
875 
876  /************************************************/
877  return;
878 }
879 
880 
881 /* >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< */
882 /* >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< */
883 /* >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< */
884 
885 /****************************************************************/
886 /* */
887 /* from followLineCase.c March 16, 2001 */
888 /* */
889 /* program to find and identify bins and borders */
890 /* in the line case */
891 /* */
892 /* Authors: Sintes, A.M */
893 /* */
894 /****************************************************************/
895 
896 
897 /****************************************************************/
898 
899 static void FollowLineCase( INT4 currentBin, UINT4 *lastBorderP,
900  REAL8 alpha, REAL8 delta, REAL8 eps, REAL8 rcOld,
901  INT4 direction, INT4 *ifailP, HOUGHptfLUT *lut,
902  HOUGHPatchGrid *patch )
903 {
904 
905  INT4 lastBorder;
906 
907 
908  REAL8 xA, yA;
909  INT4 yymin, yymax;
910  INT4 noIn; /* if no intersection occurs noIn=0 */
911 
912 
913  /* some paranoid error checking */
914  if ( patch == NULL ) {
915  fprintf( stderr, "pointer patch is null [ConstructPLUT.c %d]\n", __LINE__ );
916  }
917 
918  if ( lut == NULL ) {
919  fprintf( stderr, "pointer lut is null [ConstructPLUT.c %d]\n", __LINE__ );
920  }
921 
922  if ( ifailP == NULL ) {
923  fprintf( stderr, "pointer ifailP is null [ConstructPLUT.c %d]\n", __LINE__ );
924  }
925 
926  if ( lastBorderP == NULL ) {
927  fprintf( stderr, "pointer lastBorderP is null [ConstructPLUT.c %d]\n", __LINE__ );
928  }
929  lastBorder = *lastBorderP;
930 
931  /************************************************/
932  FindLine( alpha, delta, eps, &xA, &yA );
933  CheckLineIntersection( alpha, xA, yA, &yymin, &yymax, &noIn, patch );
934 
935  if ( noIn == 0 ) {
936  *ifailP = 0; /* =1 (ok), =0 (stop) */
937  return;
938  }
939  ++lastBorder;
940 
941  if ( yymin < 0 ) {
942  fprintf( stderr, "WARNING: Fixing yymin (%d -> 0) [ConstructPLUT.c %d]\n",
943  yymin, __LINE__ );
944  yymin = 0;
945  }
946  if ( yymax >= patch->ySide ) {
947  fprintf( stderr, "WARNING: Fixing yymax (%d -> %d) [ConstructPLUT.c %d]\n",
948  yymax, patch->ySide - 1, __LINE__ );
949  yymax = patch->ySide - 1;
950  }
951 
952  lut->border[lastBorder].yUpper = yymax;
953  lut->border[lastBorder].yLower = yymin;
954 
955  DrawLine( alpha, xA, yA, yymin, yymax,
956  &lut->border[lastBorder].xPixel[0], patch );
957 
958  /************************************************/
959  if ( direction == 1 ) {
960  /* This means that the initial case was a circle
961  that deformed into a line when increasing radius,
962  an old value of rcOld does exist */
963  REAL8 xcOld, ycOld, xRel, slope;
964 
965  xcOld = rcOld * cos( alpha );
966  ycOld = rcOld * sin( alpha );
967 
968  if ( ( alpha == LAL_PI * 0.5 ) || ( alpha == LAL_PI * 1.5 ) ||
969  ( alpha == -LAL_PI * 0.5 ) ) { /* horizontal line */
970 
971  lut->bin[currentBin].leftB1 = lastBorder;
972  lut->bin[currentBin + 1].rightB1 = lastBorder;
973 
974  /* alternatively, one can also set two extra borders */
975  /* lut->bin[currentBin].rightB2 = lastBorder; */
976  /* lut->bin[currentBin+1].leftB2 = lastBorder; */
977 
978  if ( ycOld < yA ) {
979  lut->border[lastBorder].yCenter = 0; /*or smaller*/
980  } else {
981  lut->border[lastBorder].yCenter = patch->ySide - 1; /*or bigger */
982  }
983 
984  } else { /* non- horizontal line */
985  xRel = xA + tan( alpha ) * ( yA - ycOld );
986 
987  if ( ( alpha == 0 ) || ( alpha == LAL_PI ) || ( alpha == -LAL_PI ) ) {
988  /* vertical line */
989  slope = +1; /* arbitrary number, just to avoid overflow */
990  } else {
991  slope = - cot( alpha );
992  }
993 
994  if ( xRel < xcOld ) {
995  lut->bin[currentBin].leftB1 = lastBorder;
996  lut->bin[currentBin + 1].rightB1 = lastBorder;
997 
998  if ( slope > 0 ) {
999  lut->border[lastBorder].yCenter = 0; /*or smaller*/
1000  } else {
1001  lut->border[lastBorder].yCenter = patch->ySide - 1; /*or bigger */
1002  }
1003 
1004  } else {
1005  lut->bin[currentBin].rightB2 = lastBorder;
1006  lut->bin[currentBin + 1].leftB2 = lastBorder;
1007 
1008  if ( slope < 0 ) {
1009  lut->border[lastBorder].yCenter = 0; /*or smaller*/
1010  } else {
1011  lut->border[lastBorder].yCenter = patch->ySide - 1; /*or bigger */
1012  }
1013  }
1014  } /* end non-horizontal line */
1015  } /* end of +1 direction*/
1016  else {
1017  /************************************************/
1018  /* This means ( direction == -1) */
1019  /************************************************/
1020  /* This means the first case was a line.
1021  Latter on, it can be transformed into circles,
1022  when radius decreases. No value of rcOld set */
1023 
1024  REAL8 xRel, slope;
1025 
1026  if ( ( alpha == LAL_PI * 0.5 ) || ( alpha == LAL_PI * 1.5 ) ||
1027  ( alpha == -LAL_PI * 0.5 ) ) { /* horizontal line */
1028 
1029  if ( yA < 0 ) { /* convention */
1030  lut->bin[currentBin].rightB1 = lastBorder;
1031  lut->bin[currentBin + 1].leftB1 = lastBorder;
1032  lut->border[lastBorder].yCenter = 0; /*or smaller*/
1033  } else {
1034  lut->bin[currentBin].leftB2 = lastBorder;
1035  lut->bin[currentBin + 1].rightB2 = lastBorder;
1036  lut->border[lastBorder].yCenter = patch->ySide - 1; /*or bigger */
1037  }
1038  } else { /* non- horizontal line */
1039  xRel = xA + tan( alpha ) * yA;
1040 
1041  if ( ( alpha == 0 ) || ( alpha == LAL_PI ) || ( alpha == -LAL_PI ) ) {
1042  /* vertical line */
1043  slope = +1; /* arbitrary number, just to avoid overflow */
1044  } else {
1045  slope = - cot( alpha );
1046  }
1047 
1048  if ( xRel < 0 ) {
1049  lut->bin[currentBin].leftB2 = lastBorder;
1050  lut->bin[currentBin + 1].rightB2 = lastBorder;
1051 
1052  if ( slope < 0 ) {
1053  lut->border[lastBorder].yCenter = 0; /*or smaller*/
1054  } else {
1055  lut->border[lastBorder].yCenter = patch->ySide - 1; /*or bigger */
1056  }
1057 
1058  } else {
1059  lut->bin[currentBin].rightB1 = lastBorder;
1060  lut->bin[currentBin + 1].leftB1 = lastBorder;
1061 
1062  if ( slope > 0 ) {
1063  lut->border[lastBorder].yCenter = 0; /*or smaller*/
1064  } else {
1065  lut->border[lastBorder].yCenter = patch->ySide - 1; /*or bigger */
1066  }
1067  }
1068  } /* end non-horizontal line */
1069 
1070  }/* end of -1 direction*/
1071 
1072 
1073  /************************************************/
1074 
1075  *lastBorderP = lastBorder;
1076 
1077  return;
1078 }
1079 
1080 
1081 /* >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< */
1082 /* >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< */
1083 /* >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< */
1084 
1085 /****************************************************************/
1086 /* from lines.c March 19, 2001 */
1087 /* */
1088 /* Subroutines check, find and clip lines on the patch */
1089 /* */
1090 /* */
1091 /* Authors: Sintes, A.M */
1092 /* */
1093 /****************************************************************/
1094 
1095 
1096 /* First version: to be improved thinking about
1097  substituting "ceil" and "floor" for the corresponding values
1098  assuming positive arguments. */
1099 
1100 /****************************************************************/
1101 /* Subroutine to check if we are in the line case, numerical
1102  lines inluded */
1103 /****************************************************************/
1104 
1105 static void CheckLineCase( REAL8 epsilon, REAL8 ang1, REAL8 ang2,
1106  REAL8 *epsP, INT4 *lineCaseP )
1107 {
1108 
1109  /* lineCaseP =1 line Case, =0 circle case */
1110 
1111  /* exact lines occur if */
1112  /* ((ang1==LAL_PI/2.) || (ang2==LAL_PI/2.) ||(ang2==-3.*LAL_PI/2.)) */
1113  /* other cases correspond to numerical lines */
1114 
1115  REAL8 e1, e21, e22;
1116 
1117 
1118  /* some paranoid error checking */
1119  if ( epsP == NULL ) {
1120  fprintf( stderr, "pointer epsP is null [ConstructPLUT.c %d]\n", __LINE__ );
1121  }
1122 
1123  if ( lineCaseP == NULL ) {
1124  fprintf( stderr, "pointer lineCaseP is null [ConstructPLUT.c %d]\n", __LINE__ );
1125  }
1126 
1127 
1128 
1129  e1 = ang1 - LAL_PI * 0.5;
1130  e21 = ang2 - LAL_PI * 0.5;
1131  e22 = e21 + 2 * LAL_PI;
1132 
1133  if ( fabs( e1 ) < epsilon ) {
1134  *lineCaseP = 1;
1135  *epsP = e1;
1136  return;
1137  }
1138 
1139  if ( fabs( e21 ) < epsilon ) {
1140  *lineCaseP = 1;
1141  *epsP = e21;
1142  return;
1143  }
1144 
1145  if ( fabs( e22 ) < epsilon ) {
1146  *lineCaseP = 1;
1147  *epsP = e22;
1148  return;
1149  }
1150 
1151 
1152  *lineCaseP = 0;
1153  return;
1154 }
1155 
1156 
1157 /****************************************************************/
1158 /* Subroutine to find the parameters of the exact line!!! */
1159 /* The line equation: */
1160 /* y = cotg(alpha) [xA - x ] +yA */
1161 /****************************************************************/
1162 
1163 #if 0 /* NOT USED */
1164 static void FindExactLine( REAL8 alpha, REAL8 delta,
1165  REAL8 *xA, REAL8 *yA )
1166 {
1167 
1168  REAL8 lambda, rA;
1169 
1170  lambda = 2.*delta - LAL_PI * 0.5;
1171 
1172  if ( fabs( 1. - sin( lambda ) < 1.0e-6 ) {
1173  fprintf( stderr, "warning: possible division by 0 [ConstructPLUT.c %d]\n", __LINE__ );
1174  }
1175  rA = 2.* cos( lambda ) / ( 1. - sin( lambda ) );
1176  *xA = rA * cos( alpha );
1177  *yA = rA * sin( alpha );
1178 
1179  return;
1180 }
1181 #endif
1182 
1183 /***************************************************************/
1184 /* Subroutine to find the parameters of 'approximated' lines */
1185 /* The line equation: */
1186 /* y = cotg(alpha) [xA - x ] +yA */
1187 /***************************************************************/
1188 
1189 static void FindLine( REAL8 alpha, REAL8 delta, REAL8 eps,
1190  REAL8 *xA, REAL8 *yA )
1191 {
1192 
1193  REAL8 lambda, rA;
1194 
1195  lambda = 2.*delta - LAL_PI * 0.5 - eps;
1196  rA = 2.* cos( lambda ) / ( 1. - sin( lambda ) );
1197  *xA = rA * cos( alpha );
1198  *yA = rA * sin( alpha );
1199 
1200  return;
1201 }
1202 
1203 
1204 /*********************************************************/
1205 /* Subroutine to check if the line intersects the patch. */
1206 /* Output: */
1207 /* information if intersects the patch: */
1208 /* noIn =0 no intersection */
1209 /* noIn =1 intersection */
1210 /* and the y-pixel range of the intersection */
1211 /*********************************************************/
1212 
1213 static void CheckLineIntersection( REAL8 alpha, REAL8 xA, REAL8 yA,
1214  INT4 *yyminP, INT4 *yymaxP, INT4 *noInP,
1215  HOUGHPatchGrid *patch )
1216 {
1217 
1218  INT4 yymin, yymax, noIn;
1219  volatile REAL4 kkk;
1220 
1221  yymin = 0;
1222  yymax = 0;
1223  noIn = 0;
1224 
1225  if ( ( alpha == 0 ) || ( alpha == LAL_PI ) || ( alpha == -LAL_PI ) ) {
1226  /* vertical line */
1227  if ( ( patch->xMin <= xA ) && ( patch->xMax >= xA ) ) {
1228  noIn = 1;
1229  yymin = 0;
1230  yymax = patch->ySide - 1;
1231  }
1232  } else {
1233  if ( ( alpha == LAL_PI * 0.5 ) || ( alpha == LAL_PI * 1.5 ) ||
1234  ( alpha == -LAL_PI * 0.5 ) ) {
1235  /* horizontal line */
1236  if ( ( patch->yMin <= yA ) && ( patch->yMax >= yA ) ) {
1237  noIn = 1;
1238  /* yymin = ceil((REAL4)(yA/patch->deltaY-0.5) + (REAL4)(patch->ySide*0.5)); */
1239  kkk = yA / patch->deltaY - 0.5;
1240  kkk += patch->ySide * 0.5;
1241  yymin = ceil( kkk );
1242 
1243  /* yymax = floor((REAL4)(yA/patch->deltaY-0.5) +(REAL4)(patch->ySide*0.5)); */
1244  kkk = yA / patch->deltaY - 0.5;
1245  kkk += patch->ySide * 0.5;
1246  yymax = floor( kkk );
1247  /* Note yymax < yymin, to identify an horizontal line!
1248  If the area to mark is above (arriba), mark yymin and higher.
1249  If the area to mark is below (abajo), mark yymax and lower.*/
1250  }
1251  } else {
1252  /* generic case */
1253  REAL8 myy1, y2, slope;
1254  slope = cot( alpha );
1255  myy1 = slope * ( xA - patch->xMin ) + yA;
1256  y2 = slope * ( xA - patch->xMax ) + yA;
1257 
1258  if ( ( ( myy1 >= patch->yMin ) || ( y2 >= patch->yMin ) )
1259  && ( ( myy1 <= patch->yMax ) || ( y2 <= patch->yMax ) ) ) {
1260  REAL8 yupper, ylower;
1261  noIn = 1;
1262  yupper = MAX( myy1, y2 );
1263  ylower = MIN( myy1, y2 ); /* or ylower=myy1+y2-yupper */
1264  yupper = MIN( yupper, patch->yMax );
1265  ylower = MAX( ylower, patch->yMin );
1266  /* yymin = ceil((REAL4)(ylower/patch->deltaY -0.5)+(REAL4)(patch->ySide*0.5)); */
1267  kkk = ylower / patch->deltaY - 0.5;
1268  kkk += patch->ySide * 0.5;
1269  yymin = ceil( kkk );
1270  /* yymax =floor((REAL4)(yupper/patch->deltaY-0.5)+(REAL4)(patch->ySide*0.5)); */
1271  kkk = yupper / patch->deltaY - 0.5;
1272  kkk += patch->ySide * 0.5;
1273  yymax = floor( kkk );
1274  /* in case of a pseudo-horizontal line (yymax < yymin) */
1275  /* we use the same convention as in the horizontal case */
1276  }
1277  }
1278  }
1279 
1280  *yyminP = yymin;
1281  *yymaxP = yymax;
1282  *noInP = noIn;
1283 
1284  return;
1285 }
1286 
1287 
1288 /*****************************************************************/
1289 /* Subroutine to draw non-horizontal lines! */
1290 /*****************************************************************/
1291 static void DrawLine( REAL8 alpha, REAL8 xA, REAL8 yA,
1292  INT4 yymin, INT4 yymax, COORType *column,
1293  HOUGHPatchGrid *patch )
1294 {
1295 
1296  INT4 jj;
1297  volatile REAL4 kkk;
1298 
1299  column[yymin] = patch->xSide;
1300  column[yymax] = patch->xSide;
1301 
1302  /* the if-else, is not really needed, just to avoid repeating
1303  the ceil when not necessary! */
1304 
1305  if ( ( alpha == 0 ) || ( alpha == LAL_PI ) ) {
1306  /* vertical line */
1307  INT4 xpixel;
1308 
1309  /* xpixel = ceil( (REAL4)(xA/patch->deltaX-0.5) +(REAL4)(patch->xSide*0.5)); */
1310  kkk = xA / patch->deltaX - 0.5;
1311  kkk += patch->xSide * 0.5;
1312  xpixel = ceil( kkk );
1313 
1314  if ( xpixel < 0 ) {
1315  xpixel = 0;
1316  }
1317  if ( xpixel > patch->xSide ) {
1318  xpixel = patch->xSide;
1319  }
1320 
1321  for ( jj = yymin; jj <= yymax; ++jj ) {
1322  column[jj] = xpixel;
1323  }
1324  } else {
1325  /* remaining cases */
1326  REAL8 tanalpha;
1327  REAL8 xofy;
1328  INT4 xpixel;
1329 
1330  tanalpha = tan( alpha );
1331 
1332  for ( jj = yymin; jj <= yymax; ++jj ) {
1333  /* will not be executed in the horizontal case */
1334  xofy = xA + tanalpha * ( yA - patch->yCoor[jj] );
1335  kkk = xofy / patch->deltaX - 0.5;
1336  kkk += patch->xSide * 0.5;
1337  xpixel = ceil( kkk );
1338 
1339  if ( xpixel < 0 ) {
1340  xpixel = 0;
1341  }
1342  if ( xpixel > patch->xSide ) {
1343  xpixel = patch->xSide;
1344  }
1345 
1346  column[jj] = xpixel;
1347  }
1348  }
1349 
1350  return;
1351 }
1352 
1353 
1354 
1355 /* >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< */
1356 /* >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< */
1357 /* >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< */
1358 
1359 /****************************************************************/
1360 /* */
1361 /* from initialCircleCase.c March 16, 2001 */
1362 /* */
1363 /* program to find and identify bins and borders */
1364 /* in the circle case */
1365 /* */
1366 /* Authors: Sintes, A.M */
1367 /* */
1368 /****************************************************************/
1369 
1370 
1371 
1372 /****************************************************************/
1373 
1374 static void InitialCircleCase( UINT4 *lastBorderP, REAL8 alpha,
1375  REAL8 ang1, REAL8 ang2,
1376  REAL8 *rcOldP, INT4 *directionP,
1377  INT4 *ifailP, HOUGHptfLUT *lut,
1378  HOUGHPatchGrid *patch )
1379 {
1380 
1381  INT4 lastBorder;
1382  INT4 direction; /* +1, or -1 */
1383 
1384 
1385  REAL8 rho1, rho2, radius;
1386  REAL8 xc, yc, rc; /* coordinates of the center of the circle */
1387  INT4 pieces;
1388 
1389 
1390  lastBorder = *lastBorderP;
1391 
1392  rho1 = cos( ang1 ) / ( 1. - sin( ang1 ) );
1393  rho2 = cos( ang2 ) / ( 1. - sin( ang2 ) );
1394  rc = rho1 + rho2; /* positive or negative */
1395  xc = rc * cos( alpha );
1396  yc = rc * sin( alpha );
1397  radius = fabs( rho1 - rho2 ); /* abs could be avoided */
1398 
1399  /************************************************/
1400  /* check for exclusion of intersection of the circle with the patch*/
1401  if ( ( yc + radius < patch->yMin ) || ( yc - radius > patch->yMax ) ||
1402  ( xc + radius < patch->xMin ) || ( xc - radius > patch->xMax ) ||
1403  ( sqrt( patch->xMax * patch->xMax + patch->yMax * patch->yMax ) + radius < fabs( rc ) ) ) {
1404  /* no intersection */
1405  *ifailP = 0; /* =1 (ok), =0 (stop) */
1406  return;
1407  }
1408 
1409  /************************************************/
1410  /* possible intersection case: */
1411  /************************************************/
1412 
1413  /* direction = +1 (+,-) (left, right) */
1414  /* direction = -1 (-,+) (right, left) */
1415  /************************************************/
1416  direction = ( ( fabs( rc ) < radius ) ? ( 1 ) : ( -1 ) );
1417  /************************************************/
1418 
1419 
1420  if ( direction == -1 ) {
1421  pieces = 0;
1422 
1423  /* check left circle */
1424  if ( xc >= patch->xMin ) { /* draw ( ? */
1425  INT4 yymax;
1426  INT4 yymin;
1427  INT4 noIn; /* if no intersection occurs noIn=0 */
1428 
1429  CheckLeftCircle( xc, yc, radius, &yymin, &yymax, &noIn, patch );
1430  if ( noIn ) {
1431  ++pieces;
1432  ++lastBorder;
1433 
1434  if ( yymin < 0 ) {
1435  fprintf( stderr, "WARNING: Fixing yymin (%d -> 0) [ConstructPLUT.c %d]\n",
1436  yymin, __LINE__ );
1437  yymin = 0;
1438  }
1439  if ( yymax >= patch->ySide ) {
1440  fprintf( stderr, "WARNING: Fixing yymax (%d -> %d) [ConstructPLUT.c %d]\n",
1441  yymax, patch->ySide - 1, __LINE__ );
1442  yymax = patch->ySide - 1;
1443  }
1444 
1445  lut->border[lastBorder].yUpper = yymax;
1446  lut->border[lastBorder].yLower = yymin;
1447  lut->border[lastBorder].yCenter =
1448  floor( yc / patch->deltaY + patch->ySide * 0.5 );
1449  /* rint( yc/patch->deltaY + patch->ySide*0.5 -0.5); */
1450  DrawLeftCircle( xc, yc, radius, yymin, yymax,
1451  &lut->border[lastBorder].xPixel[0], patch );
1452  lut->bin[0].rightB1 = lastBorder;
1453  lut->bin[1].leftB1 = lastBorder;
1454  }
1455  }
1456 
1457  /* check right circle */
1458  if ( xc <= patch->xMax ) { /* draw ) ? */
1459  INT4 yymax;
1460  INT4 yymin;
1461  INT4 noIn; /* if no intersection occurs noIn=0 */
1462 
1463  CheckRightCircle( xc, yc, radius, &yymin, &yymax, &noIn, patch );
1464  if ( noIn ) {
1465  ++pieces;
1466  ++lastBorder;
1467 
1468  if ( yymin < 0 ) {
1469  fprintf( stderr, "WARNING: Fixing yymin (%d -> 0) [ConstructPLUT.c %d]\n",
1470  yymin, __LINE__ );
1471  yymin = 0;
1472  }
1473  if ( yymax >= patch->ySide ) {
1474  fprintf( stderr, "WARNING: Fixing yymax (%d -> %d) [ConstructPLUT.c %d]\n",
1475  yymax, patch->ySide - 1, __LINE__ );
1476  yymax = patch->ySide - 1;
1477  }
1478 
1479  lut->border[lastBorder].yUpper = yymax;
1480  lut->border[lastBorder].yLower = yymin;
1481  lut->border[lastBorder].yCenter =
1482  floor( yc / patch->deltaY + patch->ySide * 0.5 );
1483  /*rint( yc/patch->deltaY + patch->ySide*0.5 -0.5); */
1484  DrawRightCircle( xc, yc, radius, yymin, yymax,
1485  &lut->border[lastBorder].xPixel[0], patch );
1486  lut->bin[0].leftB2 = lastBorder;
1487  lut->bin[1].rightB2 = lastBorder;
1488  }
1489  }
1490 
1491  /* real no intersection case */
1492  if ( pieces == 0 ) {
1493  *ifailP = 0;
1494  return;
1495  }
1496 
1497  } else {
1498 
1499  /************************************************/
1500  /* This means ( direction == +1) */
1501  /************************************************/
1502  /* no pathologies */
1503 
1504  pieces = 0;
1505  /* check left circle */
1506  if ( xc >= patch->xMin ) { /* draw ( ? */
1507  INT4 yymax;
1508  INT4 yymin;
1509  INT4 noIn; /* if no intersection occurs noIn=0 */
1510 
1511  CheckLeftCircle( xc, yc, radius, &yymin, &yymax, &noIn, patch );
1512  if ( noIn ) {
1513  ++pieces;
1514  ++lastBorder;
1515 
1516  if ( yymin < 0 ) {
1517  fprintf( stderr, "WARNING: Fixing yymin (%d -> 0) [ConstructPLUT.c %d]\n",
1518  yymin, __LINE__ );
1519  yymin = 0;
1520  }
1521  if ( yymax >= patch->ySide ) {
1522  fprintf( stderr, "WARNING: Fixing yymax (%d -> %d) [ConstructPLUT.c %d]\n",
1523  yymax, patch->ySide - 1, __LINE__ );
1524  yymax = patch->ySide - 1;
1525  }
1526 
1527  lut->border[lastBorder].yUpper = yymax;
1528  lut->border[lastBorder].yLower = yymin;
1529  lut->border[lastBorder].yCenter =
1530  floor( yc / patch->deltaY + patch->ySide * 0.5 );
1531  /* rint( yc/patch->deltaY + patch->ySide*0.5 -0.5); */
1532  DrawLeftCircle( xc, yc, radius, yymin, yymax,
1533  &lut->border[lastBorder].xPixel[0], patch );
1534 
1535  lut->bin[0].leftB1 = lastBorder;
1536  lut->bin[1].rightB1 = lastBorder;
1537  }
1538  }
1539 
1540  /* check right circle */
1541  if ( xc <= patch->xMax ) { /* draw ) ? */
1542  INT4 yymax;
1543  INT4 yymin;
1544  INT4 noIn; /* if no intersection occurs noIn=0 */
1545 
1546  CheckRightCircle( xc, yc, radius, &yymin, &yymax, &noIn, patch );
1547  if ( noIn ) {
1548  ++pieces;
1549  ++lastBorder;
1550 
1551  if ( yymin < 0 ) {
1552  fprintf( stderr, "WARNING: Fixing yymin (%d -> 0) [ConstructPLUT.c %d]\n",
1553  yymin, __LINE__ );
1554  yymin = 0;
1555  }
1556  if ( yymax >= patch->ySide ) {
1557  fprintf( stderr, "WARNING: Fixing yymax (%d -> %d) [ConstructPLUT.c %d]\n",
1558  yymax, patch->ySide - 1, __LINE__ );
1559  yymax = patch->ySide - 1;
1560  }
1561 
1562  lut->border[lastBorder].yUpper = yymax;
1563  lut->border[lastBorder].yLower = yymin;
1564  lut->border[lastBorder].yCenter =
1565  floor( yc / patch->deltaY + patch->ySide * 0.5 );
1566  /* rint( yc/patch->deltaY + patch->ySide*0.5 -0.5); */
1567  DrawRightCircle( xc, yc, radius, yymin, yymax,
1568  &lut->border[lastBorder].xPixel[0], patch );
1569 
1570  lut->bin[0].rightB2 = lastBorder;
1571  lut->bin[1].leftB2 = lastBorder;
1572  }
1573  }
1574 
1575  /* real no intersection case */
1576  if ( pieces == 0 ) {
1577  *ifailP = 0;
1578  return;
1579  }
1580 
1581  }
1582 
1583  /************************************************/
1584  *lastBorderP = lastBorder;
1585  *directionP = direction;
1586  *rcOldP = rc;
1587 
1588 
1589  return;
1590 }
1591 
1592 
1593 /* >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< */
1594 /* >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< */
1595 /* >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< */
1596 
1597 
1598 /****************************************************************/
1599 /* */
1600 /* from secondCircleCase.c March 16, 2001 */
1601 /* */
1602 /* program to find and identify bins and borders */
1603 /* in the circle case */
1604 /* */
1605 /* Authors: Sintes, A.M */
1606 /* */
1607 /****************************************************************/
1608 
1609 
1610 /****************************************************************/
1611 
1612 static void SecondCircleCase( INT4 currentBin, UINT4 *lastBorderP,
1613  REAL8 alpha, REAL8 ang1, REAL8 ang2,
1614  INT4 directionPlus, REAL8 *rcOldP,
1615  INT4 *pathologyP, INT4 *directionP,
1616  INT4 *ifailP, HOUGHptfLUT *lut,
1617  HOUGHPatchGrid *patch )
1618 {
1619 
1620  INT4 lastBorder;
1621  INT4 pathology; /* =1 (normal), =0 (anormal) */
1622  INT4 direction; /* +1, or -1 */
1623 
1624 
1625  REAL8 rho1, rho2, radius;
1626  REAL8 xc, yc, rc; /* coordinates of the center of the circle */
1627  INT4 pieces;
1628 
1629  pathology = *pathologyP;
1630  lastBorder = *lastBorderP;
1631 
1632  rho1 = cos( ang1 ) / ( 1. - sin( ang1 ) );
1633  rho2 = cos( ang2 ) / ( 1. - sin( ang2 ) );
1634  rc = rho1 + rho2; /* positive or negative */
1635  xc = rc * cos( alpha );
1636  yc = rc * sin( alpha );
1637  radius = fabs( rho1 - rho2 ); /* abs could be avoided */
1638 
1639  /************************************************/
1640  /* check for exclusion of intersection of the circle with the patch*/
1641  if ( ( yc + radius < patch->yMin ) || ( yc - radius > patch->yMax ) ||
1642  ( xc + radius < patch->xMin ) || ( xc - radius > patch->xMax ) ||
1643  ( sqrt( patch->xMax * patch->xMax + patch->yMax * patch->yMax ) + radius < fabs( rc ) ) ) {
1644  /* no intersection */
1645  *ifailP = 0; /* =1 (ok), =0 (stop) */
1646  return;
1647  }
1648 
1649  /************************************************/
1650  /* possible intersection case: */
1651  /************************************************/
1652 
1653  /* direction = +1 (+,-) (left, right) */
1654  /* direction = -1 (-,+) (right, left) */
1655  /************************************************/
1656  direction = ( ( fabs( rc ) < radius ) ? ( 1 ) : ( -1 ) );
1657  /************************************************/
1658 
1659 
1660  if ( direction == -1 ) { /* possible pathologies */
1661 
1662  if ( directionPlus == -1 ) {
1663  pathology = 0;
1664  }
1665 
1666  pieces = 0;
1667 
1668  /* check left circle */
1669  if ( xc >= patch->xMin ) { /* draw ( ? */
1670  INT4 yymax;
1671  INT4 yymin;
1672  INT4 noIn; /* if no intersection occurs noIn=0 */
1673 
1674  CheckLeftCircle( xc, yc, radius, &yymin, &yymax, &noIn, patch );
1675  if ( noIn ) {
1676  ++pieces;
1677  ++lastBorder;
1678 
1679  if ( yymin < 0 ) {
1680  fprintf( stderr, "WARNING: Fixing yymin (%d -> 0) [ConstructPLUT.c %d]\n",
1681  yymin, __LINE__ );
1682  yymin = 0;
1683  }
1684  if ( yymax >= patch->ySide ) {
1685  fprintf( stderr, "WARNING: Fixing yymax (%d -> %d) [ConstructPLUT.c %d]\n",
1686  yymax, patch->ySide - 1, __LINE__ );
1687  yymax = patch->ySide - 1;
1688  }
1689 
1690  lut->border[lastBorder].yUpper = yymax;
1691  lut->border[lastBorder].yLower = yymin;
1692  lut->border[lastBorder].yCenter =
1693  floor( yc / patch->deltaY + patch->ySide * 0.5 );
1694  /* rint( yc/patch->deltaY + patch->ySide*0.5 -0.5); */
1695  DrawLeftCircle( xc, yc, radius, yymin, yymax,
1696  &lut->border[lastBorder].xPixel[0], patch );
1697  if ( pathology ) {
1698  lut->bin[0].rightB1 = lastBorder;
1699  } else {
1700  lut->bin[0].rightB2 = lastBorder; /* modified */
1701  }
1702  lut->bin[currentBin + 1].leftB1 = lastBorder;
1703  }
1704  }
1705 
1706  /* check right circle */
1707  if ( xc <= patch->xMax ) { /* draw ) ? */
1708  INT4 yymax;
1709  INT4 yymin;
1710  INT4 noIn; /* if no intersection occurs noIn=0 */
1711 
1712  CheckRightCircle( xc, yc, radius, &yymin, &yymax, &noIn, patch );
1713  if ( noIn ) {
1714  ++pieces;
1715  ++lastBorder;
1716 
1717  if ( yymin < 0 ) {
1718  fprintf( stderr, "WARNING: Fixing yymin (%d -> 0) [ConstructPLUT.c %d]\n",
1719  yymin, __LINE__ );
1720  yymin = 0;
1721  }
1722  if ( yymax >= patch->ySide ) {
1723  fprintf( stderr, "WARNING: Fixing yymax (%d -> %d) [ConstructPLUT.c %d]\n",
1724  yymax, patch->ySide - 1, __LINE__ );
1725  yymax = patch->ySide - 1;
1726  }
1727 
1728  lut->border[lastBorder].yUpper = yymax;
1729  lut->border[lastBorder].yLower = yymin;
1730  lut->border[lastBorder].yCenter =
1731  floor( yc / patch->deltaY + patch->ySide * 0.5 );
1732  /* rint( yc/patch->deltaY + patch->ySide*0.5 -0.5);*/
1733  DrawRightCircle( xc, yc, radius, yymin, yymax,
1734  &lut->border[lastBorder].xPixel[0], patch );
1735  if ( pathology ) {
1736  lut->bin[0].leftB2 = lastBorder;
1737  } else {
1738  lut->bin[0].leftB1 = lastBorder; /*modified */
1739  }
1740  lut->bin[currentBin + 1].rightB2 = lastBorder;
1741  }
1742  }
1743 
1744  /* real no intersection case */
1745  if ( pieces == 0 ) {
1746  *ifailP = 0;
1747  return;
1748  }
1749 
1750  } else {
1751 
1752  /************************************************/
1753  /* This means ( direction == +1) */
1754  /************************************************/
1755  /* no pathologies */
1756 
1757  pieces = 0;
1758  /* check left circle */
1759  if ( xc >= patch->xMin ) { /* draw ( ? */
1760  INT4 yymax;
1761  INT4 yymin;
1762  INT4 noIn; /* if no intersection occurs noIn=0 */
1763 
1764  CheckLeftCircle( xc, yc, radius, &yymin, &yymax, &noIn, patch );
1765  if ( noIn ) {
1766  ++pieces;
1767  ++lastBorder;
1768 
1769  if ( yymin < 0 ) {
1770  fprintf( stderr, "WARNING: Fixing yymin (%d -> 0) [ConstructPLUT.c %d]\n",
1771  yymin, __LINE__ );
1772  yymin = 0;
1773  }
1774  if ( yymax >= patch->ySide ) {
1775  fprintf( stderr, "WARNING: Fixing yymax (%d -> %d) [ConstructPLUT.c %d]\n",
1776  yymax, patch->ySide - 1, __LINE__ );
1777  yymax = patch->ySide - 1;
1778  }
1779 
1780  lut->border[lastBorder].yUpper = yymax;
1781  lut->border[lastBorder].yLower = yymin;
1782  lut->border[lastBorder].yCenter =
1783  floor( yc / patch->deltaY + patch->ySide * 0.5 );
1784  /* rint( yc/patch->deltaY + patch->ySide*0.5 -0.5); */
1785  DrawLeftCircle( xc, yc, radius, yymin, yymax,
1786  &lut->border[lastBorder].xPixel[0], patch );
1787 
1788  lut->bin[0].leftB1 = lastBorder;
1789  lut->bin[currentBin + 1].rightB1 = lastBorder;
1790  }
1791  }
1792 
1793  /* check right circle */
1794  if ( xc <= patch->xMax ) { /* draw ) ? */
1795  INT4 yymax;
1796  INT4 yymin;
1797  INT4 noIn; /* if no intersection occurs noIn=0 */
1798 
1799  CheckRightCircle( xc, yc, radius, &yymin, &yymax, &noIn, patch );
1800  if ( noIn ) {
1801  ++pieces;
1802  ++lastBorder;
1803 
1804  if ( yymin < 0 ) {
1805  fprintf( stderr, "WARNING: Fixing yymin (%d -> 0) [ConstructPLUT.c %d]\n",
1806  yymin, __LINE__ );
1807  yymin = 0;
1808  }
1809  if ( yymax >= patch->ySide ) {
1810  fprintf( stderr, "WARNING: Fixing yymax (%d -> %d) [ConstructPLUT.c %d]\n",
1811  yymax, patch->ySide - 1, __LINE__ );
1812  yymax = patch->ySide - 1;
1813  }
1814 
1815  lut->border[lastBorder].yUpper = yymax;
1816  lut->border[lastBorder].yLower = yymin;
1817  lut->border[lastBorder].yCenter =
1818  floor( yc / patch->deltaY + patch->ySide * 0.5 );
1819  /* rint( yc/patch->deltaY + patch->ySide*0.5 -0.5);*/
1820  DrawRightCircle( xc, yc, radius, yymin, yymax,
1821  &lut->border[lastBorder].xPixel[0], patch );
1822 
1823  lut->bin[0].rightB2 = lastBorder;
1824  lut->bin[currentBin + 1].leftB2 = lastBorder;
1825  }
1826  }
1827 
1828  /* real no intersection case */
1829  if ( pieces == 0 ) {
1830  *ifailP = 0;
1831  return;
1832  }
1833 
1834  }
1835 
1836 
1837  /************************************************/
1838 
1839  *lastBorderP = lastBorder;
1840  *pathologyP = pathology;
1841  *directionP = direction;
1842  *rcOldP = rc;
1843 
1844 
1845  return;
1846 }
1847 
1848 
1849 /* >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< */
1850 /* >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< */
1851 /* >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< */
1852 
1853 
1854 /****************************************************************/
1855 /* */
1856 /* from followCircleCase.c March 16, 2001 */
1857 /* */
1858 /* program to find and identify bins and borders */
1859 /* in the circle case */
1860 /* */
1861 /* Authors: Sintes, A.M */
1862 /* */
1863 /****************************************************************/
1864 
1865 
1866 /****************************************************************/
1867 
1868 static void FollowCircleCase( INT4 currentBin, UINT4 *lastBorderP, REAL8 alpha,
1869  REAL8 ang1, REAL8 ang2, REAL8 rCritic,
1870  REAL8 rcOld, INT4 *pathologyP,
1871  INT4 *directionP, INT4 *ifailP, HOUGHptfLUT *lut,
1872  HOUGHPatchGrid *patch )
1873 {
1874 
1875  INT4 lastBorder;
1876  INT4 pathology; /* =1 (normal), =0 (anormal) */
1877  INT4 direction; /* +1, or -1 */
1878 
1879 
1880  REAL8 rho1, rho2, radius;
1881  REAL8 xc, yc, rc; /* coordinates of the center of the circle */
1882  INT4 pieces;
1883 
1884 
1885  lastBorder = *lastBorderP;
1886  pathology = *pathologyP;
1887  direction = *directionP;
1888 
1889  rho1 = cos( ang1 ) / ( 1. - sin( ang1 ) );
1890  rho2 = cos( ang2 ) / ( 1. - sin( ang2 ) );
1891  rc = rho1 + rho2; /* positive or negative */
1892  xc = rc * cos( alpha );
1893  yc = rc * sin( alpha );
1894  radius = fabs( rho1 - rho2 ); /* abs could be avoided */
1895 
1896  /************************************************/
1897  /* check for exclusion of intersection of the circle with the patch*/
1898  if ( ( yc + radius < patch->yMin ) || ( yc - radius > patch->yMax ) ||
1899  ( xc + radius < patch->xMin ) || ( xc - radius > patch->xMax ) ||
1900  ( sqrt( patch->xMax * patch->xMax + patch->yMax * patch->yMax ) + radius < fabs( rc ) ) ) {
1901  /* no intersection */
1902  *ifailP = 0; /* =1 (ok), =0 (stop) */
1903  return;
1904  }
1905 
1906  /************************************************/
1907  /* possible intersection case: */
1908  /************************************************/
1909 
1910  if ( direction == -1 ) {
1911  /* no pathological cases can happen */
1912  pieces = 0;
1913 
1914  /* check left circle */
1915  if ( xc >= patch->xMin ) { /* draw ( ? */
1916  INT4 yymax;
1917  INT4 yymin;
1918  INT4 noIn; /* if no intersection occurs noIn=0 */
1919 
1920  CheckLeftCircle( xc, yc, radius, &yymin, &yymax, &noIn, patch );
1921  if ( noIn ) {
1922  ++pieces;
1923  ++lastBorder;
1924 
1925  if ( yymin < 0 ) {
1926  fprintf( stderr, "WARNING: Fixing yymin (%d -> 0) [ConstructPLUT.c %d]\n",
1927  yymin, __LINE__ );
1928  yymin = 0;
1929  }
1930  if ( yymax >= patch->ySide ) {
1931  fprintf( stderr, "WARNING: Fixing yymax (%d -> %d) [ConstructPLUT.c %d]\n",
1932  yymax, patch->ySide - 1, __LINE__ );
1933  yymax = patch->ySide - 1;
1934  }
1935 
1936  lut->border[lastBorder].yUpper = yymax;
1937  lut->border[lastBorder].yLower = yymin;
1938  lut->border[lastBorder].yCenter =
1939  floor( yc / patch->deltaY + patch->ySide * 0.5 );
1940  /* rint( yc/patch->deltaY + patch->ySide*0.5 -0.5); */
1941  DrawLeftCircle( xc, yc, radius, yymin, yymax,
1942  &lut->border[lastBorder].xPixel[0], patch );
1943  lut->bin[currentBin].rightB1 = lastBorder;
1944  lut->bin[currentBin + 1].leftB1 = lastBorder;
1945  }
1946  }
1947 
1948  /* check right circle */
1949  if ( xc <= patch->xMax ) { /* draw ) ? */
1950  INT4 yymax;
1951  INT4 yymin;
1952  INT4 noIn; /* if no intersection occurs noIn=0 */
1953 
1954  CheckRightCircle( xc, yc, radius, &yymin, &yymax, &noIn, patch );
1955  if ( noIn ) {
1956  ++pieces;
1957  ++lastBorder;
1958 
1959  if ( yymin < 0 ) {
1960  fprintf( stderr, "WARNING: Fixing yymin (%d -> 0) [ConstructPLUT.c %d]\n",
1961  yymin, __LINE__ );
1962  yymin = 0;
1963  }
1964  if ( yymax >= patch->ySide ) {
1965  fprintf( stderr, "WARNING: Fixing yymax (%d -> %d) [ConstructPLUT.c %d]\n",
1966  yymax, patch->ySide - 1, __LINE__ );
1967  yymax = patch->ySide - 1;
1968  }
1969 
1970  lut->border[lastBorder].yUpper = yymax;
1971  lut->border[lastBorder].yLower = yymin;
1972  lut->border[lastBorder].yCenter =
1973  floor( yc / patch->deltaY + patch->ySide * 0.5 );
1974  /* rint( yc/patch->deltaY + patch->ySide*0.5 -0.5); */
1975  DrawRightCircle( xc, yc, radius, yymin, yymax,
1976  &lut->border[lastBorder].xPixel[0], patch );
1977  lut->bin[currentBin].leftB2 = lastBorder;
1978  lut->bin[currentBin + 1].rightB2 = lastBorder;
1979  }
1980  }
1981 
1982  /* real no intersection case */
1983  if ( pieces == 0 ) {
1984  *ifailP = 0;
1985  return;
1986  }
1987 
1988  } else {
1989 
1990  /************************************************/
1991  /* This means ( direction == 1) */
1992  /************************************************/
1993  /* pathologies can happen, i.e. circles can convert into lines
1994  and the two circles describing the borders of an annulus
1995  do not need to be concentric but: ( ) ( ) */
1996 
1997  /* check for pathologies */
1998  if ( ( rcOld - rCritic ) * ( rc - rCritic ) < 0 ) {
1999  pathology = 0;
2000  direction = -1; /* for the next step */
2001  }
2002 
2003  pieces = 0;
2004  /* check left circle */
2005  if ( xc >= patch->xMin ) { /* draw ( ? */
2006  INT4 yymax;
2007  INT4 yymin;
2008  INT4 noIn; /* if no intersection occurs noIn=0 */
2009 
2010  CheckLeftCircle( xc, yc, radius, &yymin, &yymax, &noIn, patch );
2011  if ( noIn ) {
2012  ++pieces;
2013  ++lastBorder;
2014 
2015  if ( yymin < 0 ) {
2016  fprintf( stderr, "WARNING: Fixing yymin (%d -> 0) [ConstructPLUT.c %d]\n",
2017  yymin, __LINE__ );
2018  yymin = 0;
2019  }
2020  if ( yymax >= patch->ySide ) {
2021  fprintf( stderr, "WARNING: Fixing yymax (%d -> %d) [ConstructPLUT.c %d]\n",
2022  yymax, patch->ySide - 1, __LINE__ );
2023  yymax = patch->ySide - 1;
2024  }
2025 
2026  lut->border[lastBorder].yUpper = yymax;
2027  lut->border[lastBorder].yLower = yymin;
2028  lut->border[lastBorder].yCenter =
2029  floor( yc / patch->deltaY + patch->ySide * 0.5 );
2030  /* rint( yc/patch->deltaY + patch->ySide*0.5 -0.5); */
2031  DrawLeftCircle( xc, yc, radius, yymin, yymax,
2032  &lut->border[lastBorder].xPixel[0], patch );
2033  if ( pathology ) {
2034  lut->bin[currentBin].leftB1 = lastBorder;
2035  lut->bin[currentBin + 1].rightB1 = lastBorder;
2036  } else {
2037  lut->bin[currentBin].rightB2 = lastBorder;
2038  lut->bin[currentBin + 1].leftB1 = lastBorder;
2039  }
2040  }
2041  }
2042 
2043  /* check right circle */
2044  if ( xc <= patch->xMax ) { /* draw ) ? */
2045  INT4 yymax;
2046  INT4 yymin;
2047  INT4 noIn; /* if no intersection occurs noIn=0 */
2048 
2049  CheckRightCircle( xc, yc, radius, &yymin, &yymax, &noIn, patch );
2050  if ( noIn ) {
2051  ++pieces;
2052  ++lastBorder;
2053 
2054  if ( yymin < 0 ) {
2055  fprintf( stderr, "WARNING: Fixing yymin (%d -> 0) [ConstructPLUT.c %d]\n",
2056  yymin, __LINE__ );
2057  yymin = 0;
2058  }
2059  if ( yymax >= patch->ySide ) {
2060  fprintf( stderr, "WARNING: Fixing yymax (%d -> %d) [ConstructPLUT.c %d]\n",
2061  yymax, patch->ySide - 1, __LINE__ );
2062  yymax = patch->ySide - 1;
2063  }
2064 
2065  lut->border[lastBorder].yUpper = yymax;
2066  lut->border[lastBorder].yLower = yymin;
2067  lut->border[lastBorder].yCenter =
2068  floor( yc / patch->deltaY + patch->ySide * 0.5 );
2069  /* rint( yc/patch->deltaY + patch->ySide*0.5 -0.5);*/
2070  DrawRightCircle( xc, yc, radius, yymin, yymax,
2071  &lut->border[lastBorder].xPixel[0], patch );
2072  if ( pathology ) {
2073  lut->bin[currentBin].rightB2 = lastBorder;
2074  lut->bin[currentBin + 1].leftB2 = lastBorder;
2075  } else {
2076  lut->bin[currentBin].leftB1 = lastBorder;
2077  lut->bin[currentBin + 1].rightB2 = lastBorder;
2078  }
2079  }
2080  }
2081 
2082  /* real no intersection case */
2083  if ( pieces == 0 ) {
2084  *ifailP = 0;
2085  return;
2086  }
2087 
2088  }
2089 
2090  /************************************************/
2091 
2092  *lastBorderP = lastBorder;
2093  *pathologyP = pathology;
2094  *directionP = direction;
2095 
2096  return;
2097 }
2098 
2099 
2100 
2101 /* >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< */
2102 /* >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< */
2103 /* >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< */
2104 
2105 
2106 /****************************************************************/
2107 /* from checkCircles.c March 30, 2001 */
2108 /* */
2109 /* subroutines to check if a circle intersects */
2110 /* the patch */
2111 /* */
2112 /* */
2113 /* Authors: Sintes, A.M */
2114 /* */
2115 /****************************************************************/
2116 
2117 
2118 /*********************************************************/
2119 /* Subroutine to check if the left part of a circle */
2120 /* intersects the patch (assuming xc>=patch->xMin) */
2121 /* Output: */
2122 /* information if intersects the patch: */
2123 /* noIn =0 no intersection */
2124 /* noIn =1 intersection */
2125 /* and the y-pixel range of the intersection */
2126 /*********************************************************/
2127 
2128 static void CheckLeftCircle( REAL8 xc, REAL8 yc, REAL8 radius,
2129  INT4 *yyminP, INT4 *yymaxP, INT4 *noInP,
2130  HOUGHPatchGrid *patch )
2131 {
2132 
2133  INT4 yymin, yymax;
2134  INT4 noIn, noIn1, noIn2;
2135  REAL8 ylower;
2136  REAL8 yupper, kkpos;
2137  volatile REAL4 kkk;
2138 
2139 
2140  /*********************************************************/
2141  noIn = 0;
2142  noIn1 = 0; /* upper quadrant */
2143  noIn2 = 0; /* lower quadrant */
2144  *noInP = noIn;
2145 
2146  if ( xc < patch->xMin ) {
2147  return; /* non optimized */
2148  }
2149 
2150  /*********************************************************/
2151  /* Reduction of the interval to look at */
2152  /*********************************************************/
2153 
2154  ylower = MAX( patch->yMin, yc - radius );
2155  yupper = MIN( patch->yMax, yc + radius );
2156 
2157  /* convert to the values of the 'near' y-pixel */
2158  /* yymax = floor((REAL4) (yupper/patch->deltaY -0.5) + (REAL4) (0.5*patch->ySide)); */
2159  kkk = yupper / patch->deltaY - 0.5;
2160  kkk += 0.5 * patch->ySide;
2161  yymax = floor( kkk );
2162 
2163  kkk = ylower / patch->deltaY - 0.5;
2164  kkk += 0.5 * patch->ySide;
2165  yymin = ceil( kkk );
2166  /* yymin = ceil( (ylower/patch->deltaY-0.5) +patch->ySide/2.); */
2167 
2168  /* NEVER try here getting the pixel center, problems looking like
2169  horizontal lines */
2170  /*
2171  * ylower = patch->yCoor[yymin];
2172  * yupper = patch->yCoor[yymax];
2173  *
2174  */
2175 
2176  if ( ( yymax < 0 ) || ( yymin > patch->ySide - 1 ) ) {
2177  return; /* non optimized */
2178  }
2179 
2180  /*********************************************************/
2181  /* looking at the upper-left quadrant */
2182  /*********************************************************/
2183 
2184  if ( yc < patch->yMax ) {
2185  REAL8 x1, myy1;
2186 
2187  /* setting the yymax value */
2188  /* in case of no over-flow problems !!! */
2189 
2190  kkpos = ( radius - ( yupper - yc ) ) * ( radius + ( yupper - yc ) );
2191  kkpos = fabs( kkpos );
2192 
2193  x1 = xc - sqrt( kkpos );
2194  /* x1 = xc - sqrt( radius*radius -(yupper-yc)*(yupper-yc)); */
2195 
2196  if ( ( x1 >= patch->xMin ) && ( x1 <= patch->xMax ) ) {
2197  noIn1 = 1; /* does intersect and yymax is the value to return! */
2198  } else {
2199  if ( x1 > patch->xMax ) {
2200  kkpos = ( radius - ( patch->xMax - xc ) ) * ( radius + ( patch->xMax - xc ) );
2201  kkpos = fabs( kkpos );
2202  myy1 = yc + sqrt( kkpos );
2203  /* myy1 = yc + sqrt( radius*radius -(patch->xMax-xc)*(patch->xMax-xc));*/
2204 
2205  if ( myy1 >= patch->yMin ) {
2206  REAL8 distance2;
2207 
2208  /* yymax = floor((REAL4) (myy1/patch->deltaY-0.5) + (REAL4) (0.5*patch->ySide)); */
2209  kkk = myy1 / patch->deltaY - 0.5;
2210  kkk += 0.5 * patch->ySide;
2211  yymax = floor( kkk );
2212  distance2 = ( xc - patch->xMax ) * ( xc - patch->xMax ) + ( yc - patch->yCoor[yymax] ) * ( yc - patch->yCoor[yymax] );
2213 
2214  if ( distance2 < radius * radius ) {
2215  noIn1 = 1;
2216  }
2217 
2218  }
2219  }
2220  }
2221  }
2222  /*********************************************************/
2223  /* looking at the lower-left quadrant */
2224  /*********************************************************/
2225 
2226  if ( yc > patch->yMin ) {
2227  REAL8 x1, myy1;
2228 
2229  /*setting the yymin value */
2230  kkpos = ( radius - ( yc - ylower ) ) * ( radius + ( yc - ylower ) );
2231  kkpos = fabs( kkpos );
2232  x1 = xc - sqrt( kkpos );
2233  /* x1 = xc - sqrt( radius*radius -(yc-ylower)*(yc-ylower)); */
2234 
2235  if ( ( x1 >= patch->xMin ) && ( x1 <= patch->xMax ) ) {
2236  noIn2 = 1; /* does intersect and yymin is the value to return! */
2237  } else {
2238  if ( x1 > patch->xMax ) {
2239  kkpos = ( radius - ( patch->xMax - xc ) ) * ( radius + ( patch->xMax - xc ) );
2240  kkpos = fabs( kkpos );
2241  myy1 = yc - sqrt( kkpos );
2242 
2243  /* myy1 = yc - sqrt( radius*radius -(patch->xMax-xc)*(patch->xMax-xc));*/
2244 
2245  if ( myy1 <= patch->yMax ) {
2246  REAL8 distance2;
2247 
2248  /* yymin = ceil((REAL4) (myy1/patch->deltaY-0.5)+(REAL4)(0.5*patch->ySide)); */
2249  kkk = myy1 / patch->deltaY - 0.5;
2250  kkk += 0.5 * patch->ySide;
2251  yymin = ceil( kkk );
2252  distance2 = ( xc - patch->xMax ) * ( xc - patch->xMax ) + ( yc - patch->yCoor[yymin] ) * ( yc - patch->yCoor[yymin] );
2253  if ( distance2 < radius * radius ) {
2254  noIn2 = 1;
2255  }
2256 
2257  }
2258  }
2259  }
2260  }
2261 
2262  /*********************************************************/
2263  /* correcting yymax, yymin values */
2264  /*********************************************************/
2265 
2266  if ( noIn1 && ( !noIn2 ) ) {
2267  REAL8 x1, myy1;
2268 
2269  noIn = 1;
2270  kkpos = ( radius - ( yc - ylower ) ) * ( radius + ( yc - ylower ) );
2271  kkpos = fabs( kkpos );
2272  x1 = xc - sqrt( kkpos );
2273  /* x1 = xc - sqrt( radius*radius -(yc-ylower)*(yc-ylower));*/
2274 
2275  if ( x1 < patch->xMin ) {
2276  /* the value yymin, needs to be set correctly */
2277  kkpos = ( radius - ( patch->xMin - xc ) ) * ( radius + ( patch->xMin - xc ) );
2278  kkpos = fabs( kkpos );
2279  myy1 = yc + sqrt( kkpos );
2280  /* myy1 = yc + sqrt( radius*radius -(patch->xMin-xc)*(patch->xMin-xc));*/
2281  /* yymin = ceil( (REAL4) (myy1/patch->deltaY-0.5)+ (REAL4)(0.5*patch->ySide) ); */
2282  kkk = myy1 / patch->deltaY - 0.5;
2283  kkk += 0.5 * patch->ySide;
2284  yymin = ceil( kkk );
2285  }
2286  }
2287 
2288  /*********************************************************/
2289 
2290  if ( ( ! noIn1 ) && noIn2 ) {
2291  /* the value yymax, needs to be set correctly */
2292  REAL8 x1, myy1;
2293 
2294  noIn = 1;
2295  kkpos = ( radius - ( yupper - yc ) ) * ( radius + ( yupper - yc ) );
2296  kkpos = fabs( kkpos );
2297  x1 = xc - sqrt( kkpos );
2298 
2299  /* x1 = xc - sqrt( radius*radius -(yupper-yc)*(yupper-yc));*/
2300  if ( x1 < patch->xMin ) {
2301  /* the value yymin, needs to be set correctly */
2302  kkpos = ( radius - ( patch->xMin - xc ) ) * ( radius + ( patch->xMin - xc ) );
2303  kkpos = fabs( kkpos );
2304  myy1 = yc - sqrt( kkpos );
2305 
2306  /*myy1 = yc - sqrt( radius*radius -(patch->xMin-xc)*(patch->xMin-xc));*/
2307  kkk = myy1 / patch->deltaY - 0.5;
2308  kkk += 0.5 * patch->ySide;
2309  yymax = floor( kkk );
2310  }
2311  }
2312 
2313  /* we could also set noIn = noIn1+noIn2; */
2314  /*********************************************************/
2315  if ( noIn1 && noIn2 ) {
2316  /* the values yymin, yymax are correct */
2317  noIn = 1;
2318  }
2319  /*********************************************************/
2320 
2321  *yyminP = yymin;
2322  *yymaxP = yymax;
2323  *noInP = noIn;
2324 
2325  return;
2326 }
2327 
2328 
2329 /*********************************************************/
2330 /* Subroutine to check if the right part of a circle */
2331 /* intersects the patch. */
2332 /* Output: */
2333 /* information if intersects the patch: */
2334 /* noIn =0 nointersection */
2335 /* noIn =1 intersection */
2336 /* and the y-pixel range of the intersection */
2337 /*********************************************************/
2338 
2339 
2340 static void CheckRightCircle( REAL8 xc, REAL8 yc, REAL8 radius,
2341  INT4 *yyminP, INT4 *yymaxP, INT4 *noInP,
2342  HOUGHPatchGrid *patch )
2343 {
2344 
2345  INT4 yymin, yymax;
2346  INT4 noIn, noIn1, noIn2;
2347  REAL8 ylower;
2348  volatile REAL4 kkk;
2349  REAL8 yupper, kkpos;
2350 
2351 
2352  /*********************************************************/
2353  noIn = 0;
2354  noIn1 = 0; /* upper quadrant */
2355  noIn2 = 0; /* lower quadrant */
2356  *noInP = noIn;
2357  if ( xc > patch->xMax ) {
2358  return; /* non optimized */
2359  }
2360 
2361  /*********************************************************/
2362  /* Reduction of the interval to look at */
2363  /*********************************************************/
2364  ylower = MAX( patch->yMin, yc - radius );
2365  yupper = MIN( patch->yMax, yc + radius );
2366 
2367  /* convert to the value of the 'near' y-pixel */
2368  /* yymax = floor((REAL4) (yupper/patch->deltaY -0.5) + (REAL4) (0.5*patch->ySide)); */
2369  kkk = yupper / patch->deltaY - 0.5;
2370  kkk += 0.5 * patch->ySide;
2371  yymax = floor( kkk );
2372 
2373  /* yymin = ceil( (ylower/patch->deltaY-0.5) +patch->ySide/2.); */
2374  kkk = ( ylower / patch->deltaY - 0.5 );
2375  kkk += 0.5 * patch->ySide;
2376  yymin = ceil( kkk );
2377 
2378  if ( ( yymax < 0 ) || ( yymin > patch->ySide - 1 ) ) {
2379  return; /* non optimized */
2380  }
2381 
2382 
2383  /*********************************************************/
2384  /* looking at the upper-right quadrant */
2385  /*********************************************************/
2386 
2387  if ( yc < patch->yMax ) {
2388  REAL8 x1, myy1;
2389 
2390  /* in case of no over-flow problems !!! */
2391 
2392  kkpos = ( radius - ( yupper - yc ) ) * ( radius + ( yupper - yc ) );
2393  kkpos = fabs( kkpos );
2394  x1 = xc + sqrt( kkpos );
2395 
2396  /* x1 = xc + sqrt( radius*radius -(yupper-yc)*(yupper-yc));*/
2397 
2398  if ( ( x1 >= patch->xMin ) && ( x1 <= patch->xMax ) ) {
2399  noIn1 = 1; /* does intersect and yymax is the value to return! */
2400  } else {
2401  if ( x1 < patch->xMin ) {
2402  kkpos = ( radius - ( patch->xMin - xc ) ) * ( radius + ( patch->xMin - xc ) );
2403  kkpos = fabs( kkpos );
2404  myy1 = yc + sqrt( kkpos );
2405 
2406  /*myy1 = yc + sqrt(radius*radius -(patch->xMin-xc)*(patch->xMin-xc));*/
2407  if ( myy1 >= patch->yMin ) {
2408  noIn1 = 1;
2409  /* yymax = floor(myy1/patch->deltaY+patch->ySide/2.-0.5);*/
2410  /* yymax = floor((REAL4) (myy1/patch->deltaY-0.5) + (REAL4) (0.5*patch->ySide)); */
2411  kkk = myy1 / patch->deltaY - 0.5;
2412  kkk += 0.5 * patch->ySide;
2413  yymax = floor( kkk );
2414  }
2415  }
2416  }
2417  }
2418 
2419  /*********************************************************/
2420  /* looking at the lower-right quadrant */
2421  /*********************************************************/
2422 
2423  if ( yc > patch->yMin ) {
2424  REAL8 x1, myy1;
2425 
2426  kkpos = ( radius - ( yc - ylower ) ) * ( radius + ( yc - ylower ) );
2427  kkpos = fabs( kkpos );
2428  x1 = xc + sqrt( kkpos );
2429 
2430  /* x1 = xc + sqrt( radius*radius -(yc-ylower)*(yc-ylower));*/
2431  if ( ( x1 >= patch->xMin ) && ( x1 <= patch->xMax ) ) {
2432  noIn2 = 1; /* does intersect and yymin is the value to return! */
2433  } else {
2434  if ( x1 < patch->xMin ) {
2435  kkpos = ( radius - ( patch->xMin - xc ) ) * ( radius + ( patch->xMin - xc ) );
2436  kkpos = fabs( kkpos );
2437  myy1 = yc - sqrt( kkpos );
2438 
2439  /* myy1 = yc - sqrt( radius*radius -(patch->xMin-xc)*(patch->xMin-xc));*/
2440  if ( myy1 <= patch->yMax ) {
2441  noIn2 = 1;
2442  /* yymin = ceil(myy1/patch->deltaY+patch->ySide/2.-0.5);*/
2443  /* yymin = ceil((REAL4) (myy1/patch->deltaY-0.5)+(REAL4)(0.5*patch->ySide)); */
2444  kkk = myy1 / patch->deltaY - 0.5;
2445  kkk += 0.5 * patch->ySide;
2446  yymin = ceil( kkk );
2447  }
2448  }
2449  }
2450  }
2451 
2452  /*********************************************************/
2453  /* correcting yymax, yymin values */
2454  /*********************************************************/
2455 
2456  if ( noIn1 && ( !noIn2 ) ) {
2457  REAL8 x1, myy1;
2458 
2459  noIn = 1;
2460  kkpos = ( radius - ( yc - ylower ) ) * ( radius + ( yc - ylower ) );
2461  kkpos = fabs( kkpos );
2462  x1 = xc + sqrt( kkpos );
2463 
2464  /* x1 = xc + sqrt( radius*radius -(yc-ylower)*(yc-ylower));*/
2465  if ( x1 > patch->xMax ) {
2466  /* the value yymin, needs to be set correctly */
2467  kkpos = ( radius - ( patch->xMax - xc ) ) * ( radius + ( patch->xMax - xc ) );
2468  kkpos = fabs( kkpos );
2469  myy1 = yc + sqrt( kkpos );
2470 
2471  /* myy1 = yc + sqrt( radius*radius -(patch->xMax-xc)*(patch->xMax-xc));*/
2472  /* yymin = ceil( myy1/patch->deltaY+patch->ySide/2.-0.5 ); */
2473  /* yymin = ceil( (REAL4) (myy1/patch->deltaY-0.5)+ (REAL4)(0.5*patch->ySide) ); */
2474  kkk = myy1 / patch->deltaY - 0.5;
2475  kkk += 0.5 * patch->ySide;
2476  yymin = ceil( kkk );
2477  }
2478  }
2479 
2480  /*********************************************************/
2481 
2482  if ( ( ! noIn1 ) && noIn2 ) {
2483  /* the value yymax, needs to be set correctly */
2484  REAL8 x1, myy1;
2485 
2486  noIn = 1;
2487  kkpos = ( radius - ( yupper - yc ) ) * ( radius + ( yupper - yc ) );
2488  kkpos = fabs( kkpos );
2489  x1 = xc + sqrt( kkpos );
2490  /* x1 = xc + sqrt( radius*radius -(yupper-yc)*(yupper-yc)); */
2491  if ( x1 > patch->xMax ) {
2492  /* the value yymin, needs to be set correctly */
2493  kkpos = ( radius - ( patch->xMax - xc ) ) * ( radius + ( patch->xMax - xc ) );
2494  kkpos = fabs( kkpos );
2495  myy1 = yc - sqrt( kkpos );
2496 
2497  /* myy1 = yc - sqrt( radius*radius -(patch->xMax-xc)*(patch->xMax-xc));*/
2498  /* yymax = floor( myy1/patch->deltaY+patch->ySide/2.-0.5 ); */
2499  /* yymax = floor( (REAL4) (myy1/patch->deltaY-0.5)+ (REAL4)(0.5*patch->ySide)); */
2500  kkk = myy1 / patch->deltaY - 0.5;
2501  kkk += 0.5 * patch->ySide;
2502  yymax = floor( kkk );
2503 
2504  }
2505  }
2506 
2507  /*********************************************************/
2508 
2509  if ( noIn1 && noIn2 ) {
2510  /* the values yymin, yymax are correct */
2511  noIn = 1;
2512  }
2513 
2514  /*********************************************************/
2515 
2516  *yyminP = yymin;
2517  *yymaxP = yymax;
2518  *noInP = noIn;
2519 
2520  return;
2521 }
2522 
2523 
2524 
2525 /* >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< */
2526 /* >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< */
2527 /* >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< */
2528 
2529 
2530 /****************************************************************/
2531 /* from drawCircles.c Marc 16, 2001 */
2532 /* */
2533 /* subroutines to clip circles on the patch */
2534 /* */
2535 /* */
2536 /* Authors: Sintes, A.M */
2537 /* */
2538 /****************************************************************/
2539 
2540 
2541 /****************************************************************/
2542 /* using sqrt!!!, */
2543 /* to be substitute by the Bresenham algorithm */
2544 /****************************************************************/
2545 
2546 static void DrawLeftCircle( REAL8 xc, REAL8 yc, REAL8 radius,
2547  INT4 yymin, INT4 yymax, COORType *column,
2548  HOUGHPatchGrid *patch )
2549 {
2550  INT4 jj;
2551 
2552  for ( jj = yymin; jj <= yymax; ++jj ) {
2553  REAL8 realx, kkpos;
2554  INT4 myindex;
2555  volatile REAL4 kkk;
2556 
2557  kkpos = ( radius - ( yc - patch->yCoor[jj] ) ) * ( radius + ( yc - patch->yCoor[jj] ) );
2558  kkpos = fabs( kkpos );
2559  realx = xc - sqrt( kkpos );
2560 
2561  /*realx = xc - sqrt( radius*radius
2562  - (yc-patch->yCoor[jj])*(yc-patch->yCoor[jj]) ); */
2563  /* myindex = ceil((REAL4) (realx/patch->deltaX-0.5)+ (REAL4)(0.5*patch->xSide) ); */
2564  kkk = realx / patch->deltaX - 0.5;
2565  kkk += 0.5 * patch->xSide;
2566  myindex = ceil( kkk );
2567  if ( myindex < 0 ) {
2568  myindex = 0;
2569  }
2570  if ( myindex > patch->xSide ) {
2571  myindex = patch->xSide;
2572  }
2573 
2574  column[jj] = myindex;
2575  }
2576 
2577  return;
2578 }
2579 
2580 
2581 /****************************************************************/
2582 /* using sqrt!!!, */
2583 /* to be substitute by the Bresenham algorithm */
2584 /****************************************************************/
2585 
2586 static void DrawRightCircle( REAL8 xc, REAL8 yc, REAL8 radius,
2587  INT4 yymin, INT4 yymax, COORType *column,
2588  HOUGHPatchGrid *patch )
2589 {
2590  INT4 jj;
2591 
2592  for ( jj = yymin; jj <= yymax; ++jj ) {
2593  REAL8 realx, kkpos;
2594  INT4 myindex;
2595  volatile REAL4 kkk;
2596 
2597  kkpos = ( radius - ( yc - patch->yCoor[jj] ) ) * ( radius + ( yc - patch->yCoor[jj] ) );
2598  kkpos = fabs( kkpos );
2599  realx = xc + sqrt( kkpos );
2600 
2601  /*realx = xc + sqrt( radius*radius
2602  - (yc-patch->yCoor[jj])*( yc-patch->yCoor[jj]) ); */
2603  /* myindex = ceil((REAL4)(realx/patch->deltaX-0.5)+(REAL4)(0.5*patch->xSide)); */
2604  kkk = realx / patch->deltaX - 0.5;
2605  kkk += 0.5 * patch->xSide;
2606  myindex = ceil( kkk );
2607  if ( myindex < 0 ) {
2608  myindex = 0;
2609  }
2610  if ( myindex > patch->xSide - 1 ) {
2611  myindex = patch->xSide;
2612  }
2613  column[jj] = myindex;
2614  }
2615 
2616  return;
2617 }
2618 
2619 
2620 /* >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< */
2621 /* >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< */
2622 /* >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< */
2623 
2624 /****************************************************************/
2625 /* from fill1Column.c March 21, 2001 */
2626 /* */
2627 /* Subroutines to correct border effects */
2628 /* */
2629 /* */
2630 /* Authors: Sintes, A.M */
2631 /* */
2632 /****************************************************************/
2633 
2634 
2635 /****************************************************************/
2636 /* Normal case: for circles and lines. NOT the pathological one */
2637 /****************************************************************/
2638 
2639 static void Fill1Column( INT4 currentBin, UINT4 *lastBorderP,
2640  HOUGHptfLUT *lut, HOUGHPatchGrid *patch )
2641 {
2642 
2643  INT4 lb1, rb1, lb2, rb2; /* The border index. If zero means that */
2644  /* it does not intersect the patch, or nothing to clip */
2645  INT4 lastBorder;
2646 
2647  lastBorder = *lastBorderP;
2648 
2649  lb1 = lut->bin[currentBin].leftB1;
2650  rb1 = lut->bin[currentBin].rightB1;
2651  lb2 = lut->bin[currentBin].leftB2;
2652  rb2 = lut->bin[currentBin].rightB2;
2653 
2654  /************************************************************/
2655  /* we want to set the values */
2656  /* lut->bin[currentBin].piece* */
2657  /* and if needed set: */
2658  /* ++lastBorder; */
2659  /* lut->bin[currentBin].rightB1 = lastBorder; */
2660  /* lut->border[lastBorder].yUpper */
2661  /* lut->border[lastBorder].yLower */
2662  /* lut->border[lastBorder].xPixel[***] = 0; */
2663  /************************************************************/
2664 
2665 
2666  /* call the proper case to correct the border effect */
2667  /* could be improved using nested if's */
2668 
2669  if ( lb1 && rb1 ) {
2670  FillCaseN1( lb1, rb1, currentBin, lut, patch );
2671  return;
2672  }
2673 
2674  if ( lb1 && !( rb1 ) && !( lb2 ) ) {
2675  FillCaseN2( lb1, currentBin, lut, patch );
2676  return;
2677  }
2678 
2679  if ( lb1 && !( rb1 ) && lb2 ) {
2680  FillCaseN3( lb1, lb2, currentBin, &lastBorder, lut, patch );
2681  *lastBorderP = lastBorder;
2682  return;
2683  }
2684 
2685  if ( !( lb1 ) && rb1 && !( rb2 ) ) {
2686  FillCaseN4( rb1, currentBin, lut, patch );
2687  return;
2688  }
2689 
2690  if ( !( lb1 ) && rb1 && rb2 ) {
2691  FillCaseN5( rb1, rb2, currentBin, lut );
2692  return;
2693  }
2694 
2695  if ( !( lb1 ) && !( rb1 ) && lb2 && rb2 ) {
2696  FillCaseN6( lb2, rb2, currentBin, lut, patch );
2697  return;
2698  }
2699 
2700  if ( !( lb1 ) && !( rb1 ) && lb2 && !( rb2 ) ) {
2701  FillCaseN7( lb2, currentBin, lut, patch );
2702  return;
2703  }
2704 
2705  if ( !( lb1 ) && !( rb1 ) && !( lb2 ) && rb2 ) {
2706  FillCaseN8( rb2, currentBin, lut, patch );
2707  return;
2708  }
2709 
2710  if ( !( currentBin ) && !( lb1 ) && !( rb1 ) && !( lb2 ) && !( rb2 ) ) {
2711  lut->bin[0].piece1max = patch->ySide - 1;
2712  return;
2713  }
2714 
2715 
2716  return;
2717 }
2718 
2719 
2720 /***************************************************************/
2721 /* case: ( lb1 && rb1 ) */
2722 /***************************************************************/
2723 static void FillCaseN1( INT4 lb1, INT4 rb1, INT4 currentBin,
2724  HOUGHptfLUT *lut, HOUGHPatchGrid *patch )
2725 {
2726  INT4 lb1UpY, lb1LoY;
2727  INT4 rb1UpY, rb1LoY;
2728  INT4 yCl, yCr;
2729 
2730  lb1UpY = lut->border[lb1].yUpper;
2731  lb1LoY = lut->border[lb1].yLower;
2732 
2733  if ( ( lb1UpY == patch->ySide - 1 ) && ( lb1LoY == 0 ) ) {
2734  return;
2735  }
2736 
2737  rb1UpY = lut->border[rb1].yUpper;
2738  rb1LoY = lut->border[rb1].yLower;
2739 
2740  yCl = lut->border[lb1].yCenter;
2741  yCr = lut->border[rb1].yCenter;
2742 
2743  if ( lb1LoY > yCl ) {
2744  lut->bin[currentBin].piece1max = lb1LoY - 1;
2745  if ( rb1LoY > yCr ) {
2746  lut->bin[currentBin].piece1min = rb1LoY;
2747  }
2748  /* else{ */ /* already initialized */
2749  /* lut->bin[currentBin].piece1min = 0; } */
2750  return;
2751  }
2752 
2753  if ( lb1UpY < yCl ) {
2754  lut->bin[currentBin].piece1min = lb1UpY + 1;
2755  if ( rb1UpY < yCr ) {
2756  lut->bin[currentBin].piece1max = rb1UpY;
2757  } else {
2758  lut->bin[currentBin].piece1max = patch->ySide - 1;
2759  }
2760  }
2761 
2762  return;
2763 }
2764 
2765 /***************************************************************/
2766 /* case: (lb1 && !(rb1) && !(lb2) ) */
2767 /***************************************************************/
2768 static void FillCaseN2( INT4 lb1, INT4 currentBin,
2769  HOUGHptfLUT *lut, HOUGHPatchGrid *patch )
2770 {
2771  INT4 lb1UpY, lb1LoY;
2772  INT4 yCl;
2773 
2774  lb1UpY = lut->border[lb1].yUpper;
2775  lb1LoY = lut->border[lb1].yLower;
2776  yCl = lut->border[lb1].yCenter;
2777 
2778  if ( lb1LoY > yCl ) {
2779  lut->bin[currentBin].piece1max = lb1LoY - 1;
2780  /* lut->bin[currentBin].piece1min = 0; */
2781  return;
2782  }
2783 
2784  if ( lb1UpY < yCl ) {
2785  lut->bin[currentBin].piece1max = patch->ySide - 1;
2786  lut->bin[currentBin].piece1min = lb1UpY + 1;
2787  }
2788  return;
2789 }
2790 
2791 /***************************************************************/
2792 /* case: (lb1 && !(rb1) && lb2 ) */
2793 /***************************************************************/
2794 static void FillCaseN3( INT4 lb1, INT4 lb2, INT4 currentBin, INT4 *lastBorderP,
2795  HOUGHptfLUT *lut, HOUGHPatchGrid *patch )
2796 {
2797  INT4 lastBorder;
2798  INT4 lb1UpY, lb1LoY;
2799  INT4 lb2UpY, lb2LoY;
2800  INT4 yCl;
2801  INT4 k;
2802 
2803  lastBorder = *lastBorderP;
2804  lb1UpY = lut->border[lb1].yUpper;
2805  lb1LoY = lut->border[lb1].yLower;
2806  yCl = lut->border[lb1].yCenter;
2807  lb2UpY = lut->border[lb2].yUpper;
2808  lb2LoY = lut->border[lb2].yLower;
2809 
2810  if ( lb1LoY > yCl ) {
2811  lut->bin[currentBin].piece1max = lb1LoY - 1;
2812  lut->bin[currentBin].piece1min = lb2UpY + 1;
2813  lut->bin[currentBin].piece2max = lb2LoY - 1;
2814  /* lut->bin[currentBin].piece2min = 0; */
2815  return;
2816  }
2817 
2818  if ( lb1UpY < yCl ) {
2819  lut->bin[currentBin].piece1max = patch->ySide - 1;
2820  lut->bin[currentBin].piece1min = lb2UpY + 1;
2821  lut->bin[currentBin].piece2max = lb2LoY - 1;
2822  lut->bin[currentBin].piece2min = lb1UpY + 1;
2823  return;
2824  }
2825 
2826  /* adding a ficticious border rb1 to correct the lb1 effects */
2827  ++lastBorder;
2828  lut->bin[currentBin].rightB1 = lastBorder;
2829  lut->border[lastBorder].yUpper = lb2UpY;
2830  lut->border[lastBorder].yLower = lb2LoY;
2831  for ( k = lb2LoY; k <= lb2UpY; ++k ) {
2832  lut->border[lastBorder].xPixel[k] = 0;
2833  }
2834  *lastBorderP = lastBorder;
2835 
2836  return;
2837 }
2838 
2839 /***************************************************************/
2840 /* case: (!(lb1) && rb1 && !(rb2) ) */
2841 /***************************************************************/
2842 static void FillCaseN4( INT4 rb1, INT4 currentBin,
2843  HOUGHptfLUT *lut, HOUGHPatchGrid *patch )
2844 {
2845  INT4 rb1UpY, rb1LoY;
2846  INT4 yCr;
2847 
2848  rb1UpY = lut->border[rb1].yUpper;
2849  rb1LoY = lut->border[rb1].yLower;
2850  yCr = lut->border[rb1].yCenter;
2851 
2852  if ( rb1UpY < yCr ) {
2853  lut->bin[currentBin].piece1max = rb1UpY;
2854  /* lut->bin[currentBin].piece1min = 0; */
2855  return;
2856  }
2857 
2858  if ( rb1LoY > yCr ) {
2859  lut->bin[currentBin].piece1max = patch->ySide - 1;
2860  lut->bin[currentBin].piece1min = rb1LoY;
2861  return;
2862  }
2863 
2864  lut->bin[currentBin].piece1max = patch->ySide - 1;
2865  /* lut->bin[currentBin].piece1min = 0; */
2866 
2867  return;
2868 }
2869 
2870 /***************************************************************/
2871 /* case: (!(lb1) && rb1 && rb2 ) */
2872 /***************************************************************/
2873 static void FillCaseN5( INT4 rb1, INT4 rb2, INT4 currentBin,
2874  HOUGHptfLUT *lut )
2875 {
2876  INT4 rb1UpY, rb1LoY;
2877  INT4 rb2UpY, rb2LoY;
2878  INT4 yCr;
2879 
2880  rb1UpY = lut->border[rb1].yUpper;
2881  rb1LoY = lut->border[rb1].yLower;
2882  yCr = lut->border[rb1].yCenter;
2883  rb2UpY = lut->border[rb2].yUpper;
2884  rb2LoY = lut->border[rb2].yLower;
2885 
2886  if ( rb1UpY < yCr ) {
2887  lut->bin[currentBin].piece1max = rb1UpY;
2888  lut->bin[currentBin].piece1min = rb2LoY;
2889  return;
2890  }
2891 
2892  if ( rb1LoY > yCr ) {
2893  lut->bin[currentBin].piece1max = rb2UpY;
2894  lut->bin[currentBin].piece1min = rb1LoY;
2895  return;
2896  }
2897 
2898  lut->bin[currentBin].piece1max = rb2UpY;
2899  lut->bin[currentBin].piece1min = rb2LoY;
2900 
2901  return;
2902 }
2903 
2904 
2905 /***************************************************************/
2906 /* case: (!(lb1) && !(rb1) && lb2 && rb2) */
2907 /***************************************************************/
2908 static void FillCaseN6( INT4 lb2, INT4 rb2, INT4 currentBin,
2909  HOUGHptfLUT *lut, HOUGHPatchGrid *patch )
2910 {
2911  INT4 lb2UpY, lb2LoY;
2912  INT4 rb2UpY, rb2LoY;
2913  INT4 yCl, yCr;
2914 
2915  lb2UpY = lut->border[lb2].yUpper;
2916  lb2LoY = lut->border[lb2].yLower;
2917 
2918  if ( ( lb2UpY == patch->ySide - 1 ) && ( lb2LoY == 0 ) ) {
2919  return;
2920  }
2921 
2922  rb2UpY = lut->border[rb2].yUpper;
2923  rb2LoY = lut->border[rb2].yLower;
2924  yCl = lut->border[lb2].yCenter;
2925  yCr = lut->border[rb2].yCenter;
2926 
2927  if ( lb2UpY >= yCl ) {
2928  lut->bin[currentBin].piece1min = lb2UpY + 1;
2929  if ( rb2UpY >= yCr ) {
2930  lut->bin[currentBin].piece1max = rb2UpY;
2931  } else {
2932  lut->bin[currentBin].piece1max = patch->ySide - 1;
2933  }
2934  }
2935 
2936  if ( lb2LoY <= yCl ) {
2937  lut->bin[currentBin].piece2max = lb2LoY - 1;
2938  if ( rb2LoY <= yCr ) {
2939  lut->bin[currentBin].piece2min = rb2LoY;
2940  }
2941  /* else{ */
2942  /* lut->bin[currentBin].piece2min = 0; } */
2943  }
2944 
2945  return;
2946 }
2947 
2948 /***************************************************************/
2949 /* case: (!(lb1) && !(rb1) && lb2 && !(rb2)) */
2950 /***************************************************************/
2951 static void FillCaseN7( INT4 lb2, INT4 currentBin,
2952  HOUGHptfLUT *lut, HOUGHPatchGrid *patch )
2953 {
2954  INT4 lb2UpY, lb2LoY;
2955  INT4 yCl;
2956 
2957  lb2UpY = lut->border[lb2].yUpper;
2958  lb2LoY = lut->border[lb2].yLower;
2959  yCl = lut->border[lb2].yCenter;
2960 
2961  if ( lb2UpY >= yCl ) {
2962  lut->bin[currentBin].piece1max = patch->ySide - 1;
2963  lut->bin[currentBin].piece1min = lb2UpY + 1;
2964  }
2965 
2966  if ( lb2LoY <= yCl ) {
2967  lut->bin[currentBin].piece2max = lb2LoY - 1;
2968  /* lut->bin[currentBin].piece2min = 0; */
2969  }
2970 
2971  return;
2972 }
2973 
2974 /***************************************************************/
2975 /* case: (!(lb1) && !(rb1) && !(lb2) && rb2) */
2976 /***************************************************************/
2977 static void FillCaseN8( INT4 rb2, INT4 currentBin,
2978  HOUGHptfLUT *lut, HOUGHPatchGrid *patch )
2979 {
2980  INT4 rb2UpY, rb2LoY;
2981  INT4 yCr;
2982 
2983  rb2UpY = lut->border[rb2].yUpper;
2984  rb2LoY = lut->border[rb2].yLower;
2985  yCr = lut->border[rb2].yCenter;
2986 
2987  lut->bin[currentBin].piece1max = patch->ySide - 1;
2988  /* lut->bin[currentBin].piece1min = 0; */
2989 
2990  if ( rb2UpY >= yCr ) {
2991  lut->bin[currentBin].piece1max = rb2UpY;
2992  }
2993 
2994  if ( rb2LoY <= yCr ) {
2995  lut->bin[currentBin].piece1min = rb2LoY;
2996  }
2997 
2998  return;
2999 }
3000 
3001 
3002 
3003 /* >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< */
3004 /* >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< */
3005 /* >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< */
3006 
3007 
3008 /****************************************************************/
3009 /* from fill1ColumnAnormal.c March 23, 2001 */
3010 /* */
3011 /* Subroutines to correct border effects */
3012 /* */
3013 /* */
3014 /* Authors: Sintes, A.M */
3015 /* */
3016 /****************************************************************/
3017 
3018 
3019 /****************************************************************/
3020 /* Anormal case: for circles and lines. The pathological case */
3021 /****************************************************************/
3022 
3023 static void Fill1ColumnAnor( INT4 currentBin, HOUGHptfLUT *lut,
3024  HOUGHPatchGrid *patch )
3025 {
3026 
3027  /* It will be similar to the non-pathological case.
3028  Note that here the circles have very large radius.
3029  The y.coor of the center is very large or set to the border.
3030  Logic not defined yet! */
3031 
3032 
3033  INT4 lb1, rb1, lb2, rb2; /* The border index. If zero means that */
3034  /* it does not intersect the patch, or nothing to clip */
3035 
3036  /* The pathological case implies the name convention:
3037  lb1 is equivalent to lb2 ), rb2 equivalent to rb1 (,
3038  just to avoid border having the same name */
3039 
3040 
3041  lb1 = lut->bin[currentBin].leftB1;
3042  rb1 = lut->bin[currentBin].rightB1;
3043  lb2 = lut->bin[currentBin].leftB2;
3044  rb2 = lut->bin[currentBin].rightB2;
3045 
3046  /************************************************************/
3047  /* we want to set the values */
3048  /* lut->bin[currentBin].piece* */
3049  /************************************************************/
3050 
3051 
3052 
3053  /* call the proper case to correct the border effect */
3054 
3055  if ( rb1 && !( rb2 ) && !( lb1 ) ) {
3056  FillCaseN4( rb1, currentBin, lut, patch );
3057  return;
3058  }
3059 
3060  if ( !( rb1 ) && lb2 && !( rb2 ) && !( lb1 ) ) {
3061  FillCaseN7( lb2, currentBin, lut, patch );
3062  return;
3063  }
3064 
3065  if ( !( rb1 ) && !( lb2 ) && rb2 ) {
3066  FillCaseN4( rb2, currentBin, lut, patch );
3067  return;
3068  }
3069 
3070  if ( !( rb1 ) && !( lb2 ) && !( rb2 ) && lb1 ) {
3071  FillCaseN7( lb1, currentBin, lut, patch );
3072  return;
3073  }
3074  /************************************************************/
3075 
3076  if ( rb1 && rb2 ) {
3077  FillCaseA1( rb1, rb2, currentBin, lut );
3078  return;
3079  }
3080 
3081  if ( !( rb1 ) && lb2 && !( rb2 ) && lb1 ) {
3082  FillCaseA2( lb1, lb2, currentBin, lut );
3083  return;
3084  }
3085 
3086  if ( !( rb2 ) && lb1 && rb1 ) {
3087  FillCaseA3( lb1, rb1, currentBin, lut, patch );
3088  return;
3089  }
3090 
3091  if ( !( rb1 ) && lb2 && rb2 ) {
3092  FillCaseA3( lb2, rb2, currentBin, lut, patch );
3093  return;
3094  }
3095 
3096  if ( !( currentBin ) && !( lb1 ) && !( rb1 ) && !( lb2 ) && !( rb2 ) ) {
3097  lut->bin[0].piece1max = patch->ySide - 1;
3098  return;
3099  }
3100 
3101  return;
3102 }
3103 
3104 /***************************************************************/
3105 /* case: ( rb1 && rb2 ) */
3106 /***************************************************************/
3107 static void FillCaseA1( INT4 rb1, INT4 rb2, INT4 currentBin,
3108  HOUGHptfLUT *lut )
3109 {
3110  INT4 rb1UpY, rb1LoY;
3111  INT4 rb2UpY, rb2LoY;
3112  INT4 yCr1;
3113 
3114  rb1UpY = lut->border[rb1].yUpper;
3115  rb1LoY = lut->border[rb1].yLower;
3116  yCr1 = lut->border[rb1].yCenter;
3117 
3118  rb2UpY = lut->border[rb2].yUpper;
3119  rb2LoY = lut->border[rb2].yLower;
3120 
3121  /* since the circles are big, number of case are simplifyed */
3122  /* more cases will appear in a general case with small radius */
3123 
3124  if ( yCr1 < rb1UpY ) {
3125  lut->bin[currentBin].piece1max = rb2UpY;
3126  lut->bin[currentBin].piece1min = rb1LoY;
3127  } else {
3128  lut->bin[currentBin].piece1max = rb1UpY;
3129  lut->bin[currentBin].piece1min = rb2LoY;
3130  }
3131 
3132  return;
3133 }
3134 
3135 /***************************************************************/
3136 /* case: ( !(rb1)&& lb2 && !(rb2) && lb1 ) */
3137 /***************************************************************/
3138 static void FillCaseA2( INT4 lb1, INT4 lb2, INT4 currentBin,
3139  HOUGHptfLUT *lut )
3140 {
3141  INT4 lb1UpY, lb1LoY;
3142  INT4 lb2UpY, lb2LoY;
3143  INT4 yCl1;
3144 
3145  lb1UpY = lut->border[lb1].yUpper;
3146  lb1LoY = lut->border[lb1].yLower;
3147  yCl1 = lut->border[lb1].yCenter;
3148 
3149  lb2UpY = lut->border[lb2].yUpper;
3150  lb2LoY = lut->border[lb2].yLower;
3151 
3152  /* since the circles are big, number of case are simplifyed */
3153  /* more cases will appear in the general case */
3154 
3155  if ( yCl1 < lb1UpY ) {
3156  lut->bin[currentBin].piece1max = lb2LoY - 1;
3157  lut->bin[currentBin].piece1min = lb1UpY + 1;
3158  } else {
3159  lut->bin[currentBin].piece1max = lb1LoY - 1;
3160  lut->bin[currentBin].piece1min = lb2UpY + 1;
3161  }
3162 
3163  return;
3164 }
3165 
3166 /***************************************************************/
3167 /* case: ( !(rb2) && lb1 && rb1 ) or ( !(rb1)&& lb2 && rb2 ) */
3168 /***************************************************************/
3169 static void FillCaseA3( INT4 lb1, INT4 rb1, INT4 currentBin,
3170  HOUGHptfLUT *lut, HOUGHPatchGrid *patch )
3171 {
3172 
3173  /* here we should code all possible cases with no exceptions */
3174  INT4 rb1UpY, rb1LoY;
3175  INT4 lb1UpY, lb1LoY;
3176  INT4 yCr1, yCl1;
3177 
3178  rb1UpY = lut->border[rb1].yUpper;
3179  rb1LoY = lut->border[rb1].yLower;
3180  yCr1 = lut->border[rb1].yCenter;
3181 
3182  lb1UpY = lut->border[lb1].yUpper;
3183  lb1LoY = lut->border[lb1].yLower;
3184  yCl1 = lut->border[lb1].yCenter;
3185 
3186 
3187  if ( ( lb1UpY == patch->ySide - 1 ) && ( lb1LoY == 0 ) ) {
3188  return;
3189  }
3190 
3191  /* provisional settings */
3192  if ( lb1UpY >= yCl1 ) {
3193  lut->bin[currentBin].piece1max = patch->ySide - 1;
3194  lut->bin[currentBin].piece1min = lb1UpY + 1;
3195  }
3196 
3197  if ( lb1LoY <= yCl1 ) {
3198  lut->bin[currentBin].piece2max = lb1LoY - 1;
3199  /* lut->bin[currentBin].piece2min = 0; */
3200  }
3201 
3202  /* corrections */
3203  if ( rb1LoY > yCr1 ) {
3204  lut->bin[currentBin].piece2min = rb1LoY;
3205  return;
3206  }
3207 
3208  if ( rb1UpY < yCr1 ) {
3209  lut->bin[currentBin].piece1max = rb1UpY;
3210  }
3211 
3212  return;
3213 
3214 }
3215 
3216 /* >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< */
3217 /* >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< */
3218 /* >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< */
3219 
3220 
3221 
3222 
3223 #undef MIN
3224 #undef MAX
eps
const REAL8 eps
Definition: CWMakeFakeData.c:56
FollowLineCase
static void FollowLineCase(INT4, UINT4 *, REAL8, REAL8, REAL8, REAL8, INT4, INT4 *, HOUGHptfLUT *, HOUGHPatchGrid *)
Definition: ConstructPLUT.c:899
InitialLineCase
static void InitialLineCase(UINT4 *, REAL8, REAL8, REAL8, INT4 *, HOUGHptfLUT *, HOUGHPatchGrid *)
Definition: ConstructPLUT.c:626
FillCaseN7
static void FillCaseN7(INT4, INT4, HOUGHptfLUT *, HOUGHPatchGrid *)
Definition: ConstructPLUT.c:2951
FillCaseA2
static void FillCaseA2(INT4, INT4, INT4, HOUGHptfLUT *)
Definition: ConstructPLUT.c:3138
CheckLeftCircle
static void CheckLeftCircle(REAL8, REAL8, REAL8, INT4 *, INT4 *, INT4 *, HOUGHPatchGrid *)
Definition: ConstructPLUT.c:2128
Fill1ColumnAnor
static void Fill1ColumnAnor(INT4, HOUGHptfLUT *, HOUGHPatchGrid *)
Definition: ConstructPLUT.c:3023
Fill1Column
static void Fill1Column(INT4, UINT4 *, HOUGHptfLUT *, HOUGHPatchGrid *)
Definition: ConstructPLUT.c:2639
FillCaseN4
static void FillCaseN4(INT4, INT4, HOUGHptfLUT *, HOUGHPatchGrid *)
Definition: ConstructPLUT.c:2842
FillCaseN6
static void FillCaseN6(INT4, INT4, INT4, HOUGHptfLUT *, HOUGHPatchGrid *)
Definition: ConstructPLUT.c:2908
FollowCircleCase
static void FollowCircleCase(INT4, UINT4 *, REAL8, REAL8, REAL8, REAL8, REAL8, INT4 *, INT4 *, INT4 *, HOUGHptfLUT *, HOUGHPatchGrid *)
Definition: ConstructPLUT.c:1868
PLUTInitialize
static void PLUTInitialize(HOUGHptfLUT *)
Definition: ConstructPLUT.c:244
FillCaseN2
static void FillCaseN2(INT4, INT4, HOUGHptfLUT *, HOUGHPatchGrid *)
Definition: ConstructPLUT.c:2768
SecondCircleCase
static void SecondCircleCase(INT4, UINT4 *, REAL8, REAL8, REAL8, INT4, REAL8 *, INT4 *, INT4 *, INT4 *, HOUGHptfLUT *, HOUGHPatchGrid *)
Definition: ConstructPLUT.c:1612
DrawLine
static void DrawLine(REAL8, REAL8, REAL8, INT4, INT4, COORType *, HOUGHPatchGrid *)
Definition: ConstructPLUT.c:1291
FillCaseN5
static void FillCaseN5(INT4, INT4, INT4, HOUGHptfLUT *)
Definition: ConstructPLUT.c:2873
DrawRightCircle
static void DrawRightCircle(REAL8, REAL8, REAL8, INT4, INT4, COORType *, HOUGHPatchGrid *)
Definition: ConstructPLUT.c:2586
CheckLineCase
static void CheckLineCase(REAL8, REAL8, REAL8, REAL8 *, INT4 *)
Definition: ConstructPLUT.c:1105
FillCaseN3
static void FillCaseN3(INT4, INT4, INT4, INT4 *, HOUGHptfLUT *, HOUGHPatchGrid *)
Definition: ConstructPLUT.c:2794
FillCaseN1
static void FillCaseN1(INT4, INT4, INT4, HOUGHptfLUT *, HOUGHPatchGrid *)
Definition: ConstructPLUT.c:2723
FillCaseN8
static void FillCaseN8(INT4, INT4, HOUGHptfLUT *, HOUGHPatchGrid *)
Definition: ConstructPLUT.c:2977
InitialCircleCase
static void InitialCircleCase(UINT4 *, REAL8, REAL8, REAL8, REAL8 *, INT4 *, INT4 *, HOUGHptfLUT *, HOUGHPatchGrid *)
Definition: ConstructPLUT.c:1374
DrawLeftCircle
static void DrawLeftCircle(REAL8, REAL8, REAL8, INT4, INT4, COORType *, HOUGHPatchGrid *)
Definition: ConstructPLUT.c:2546
SecondLineCase
static void SecondLineCase(INT4, UINT4 *, REAL8, REAL8, REAL8, INT4 *, HOUGHptfLUT *, HOUGHPatchGrid *)
Definition: ConstructPLUT.c:759
FindLine
static void FindLine(REAL8, REAL8, REAL8, REAL8 *, REAL8 *)
Definition: ConstructPLUT.c:1189
FillCaseA1
static void FillCaseA1(INT4, INT4, INT4, HOUGHptfLUT *)
Definition: ConstructPLUT.c:3107
CheckLineIntersection
static void CheckLineIntersection(REAL8, REAL8, REAL8, INT4 *, INT4 *, INT4 *, HOUGHPatchGrid *)
Definition: ConstructPLUT.c:1213
FillCaseA3
static void FillCaseA3(INT4, INT4, INT4, HOUGHptfLUT *, HOUGHPatchGrid *)
Definition: ConstructPLUT.c:3169
FillPLUT
static void FillPLUT(HOUGHParamPLUT *, HOUGHptfLUT *, HOUGHPatchGrid *)
Definition: ConstructPLUT.c:332
CheckRightCircle
static void CheckRightCircle(REAL8, REAL8, REAL8, INT4 *, INT4 *, INT4 *, HOUGHPatchGrid *)
Definition: ConstructPLUT.c:2340
MAX
#define MAX(x, y)
Definition: DopplerFullScan.c:43
k
int k
delta
static double double delta
ABORT
#define ABORT(statusptr, code, mesg)
INITSTATUS
#define INITSTATUS(statusptr)
RETURN
#define RETURN(statusptr)
fprintf
#define fprintf
MIN
#define MIN(x, y)
Definition: dumpSFT.c:39
rho2
const double rho2
LAL_PI
#define LAL_PI
REAL8
double REAL8
INT8
int64_t INT8
UINT4
uint32_t UINT4
INT4
int32_t INT4
REAL4
float REAL4
LUTH_MSGENULL
#define LUTH_MSGENULL
Definition: LUT.h:157
LUTH_EVAL
#define LUTH_EVAL
Definition: LUT.h:155
LUTH_MSGESIZE
#define LUTH_MSGESIZE
Definition: LUT.h:158
cot
#define cot(A)
Definition: LUT.h:175
LUTH_ENULL
#define LUTH_ENULL
Definition: LUT.h:149
LUTH_MSGEVAL
#define LUTH_MSGEVAL
Definition: LUT.h:163
LUTH_ESIZE
#define LUTH_ESIZE
Definition: LUT.h:150
LALHOUGHConstructPLUT
void LALHOUGHConstructPLUT(LALStatus *status, HOUGHptfLUT *lut, HOUGHPatchGrid *patch, HOUGHParamPLUT *par)
Definition: ConstructPLUT.c:149
COORType
INT2 COORType
To be changed to {INT2 COORType} if the number of pixels in the x-direction exceeds 255.
Definition: LUT.h:218
lalpulsar.git_version.status
string status
Definition: git_version.py:32
lalpulsar_knope_result_page.par
par
Definition: lalpulsar_knope_result_page.py:2439
make_frame_cache.i
int i
Definition: make_frame_cache.py:77
alpha
double alpha
HOUGHBin2Border::rightB2
INT2 rightB2
Border index to be used (stop-border ‘ ’)
Definition: LUT.h:238
HOUGHBin2Border::leftB1
INT2 leftB1
Border index to be used (start-border ‘ ’)
Definition: LUT.h:235
HOUGHBin2Border::piece1min
INT2 piece1min
If piece1min piece1max no corrections should be added.
Definition: LUT.h:240
HOUGHBin2Border::piece2max
INT2 piece2max
Interval limits of the (second piece) correction to the first column.
Definition: LUT.h:241
HOUGHBin2Border::rightB1
INT2 rightB1
Border index to be used (stop-border ‘ ’)
Definition: LUT.h:236
HOUGHBin2Border::piece1max
INT2 piece1max
Interval limits of the (first piece) correction to the first column.
Definition: LUT.h:239
HOUGHBin2Border::leftB2
INT2 leftB2
Border index to be used (start-border ‘ ’)
Definition: LUT.h:237
HOUGHBin2Border::piece2min
INT2 piece2min
If piece2min piece2max no corrections should be added.
Definition: LUT.h:242
HOUGHBorder::xPixel
COORType * xPixel
x pixel index to be marked
Definition: LUT.h:226
HOUGHBorder::yLower
INT4 yLower
lower y pixel affected by this border and yUpper<yLower or yUpper<0 are possible
Definition: LUT.h:223
HOUGHBorder::yCenter
INT4 yCenter
y pixel value of the center of the circle
Definition: LUT.h:224
HOUGHBorder::yUpper
INT4 yUpper
upper y pixel affected by this border
Definition: LUT.h:222
HOUGHParamPLUT
Parameters needed to construct the partial look up table.
Definition: LUT.h:333
HOUGHPatchGrid
This structure stores patch-frequency grid information.
Definition: LUT.h:264
HOUGHPatchGrid::ySide
UINT2 ySide
Real number of pixels in the y-direction (in the projected plane).
Definition: LUT.h:275
HOUGHPatchGrid::deltaY
REAL8 deltaY
Longitudinal space resolution, y-direction.
Definition: LUT.h:272
HOUGHPatchGrid::yMin
REAL8 yMin
Patch limit, as center of the first pixel.
Definition: LUT.h:273
HOUGHPatchGrid::xSide
UINT2 xSide
Real number of pixels in the x direction (in the projected plane); it should be less than or equal to...
Definition: LUT.h:270
HOUGHPatchGrid::yMax
REAL8 yMax
Patch limit, as center of the last pixel.
Definition: LUT.h:274
HOUGHPatchGrid::deltaX
REAL8 deltaX
Longitudinal space resolution, x-direction.
Definition: LUT.h:267
HOUGHPatchGrid::xMin
REAL8 xMin
Patch limit, as the coordinate of the center of the first pixel.
Definition: LUT.h:268
HOUGHPatchGrid::xMax
REAL8 xMax
Patch limit, as the coordinate of the center of the last pixel.
Definition: LUT.h:269
HOUGHPatchGrid::yCoor
REAL8 * yCoor
Coordinates of the pixel centers.
Definition: LUT.h:278
HOUGHPatchGrid::deltaF
REAL8 deltaF
Frequency resolution: df=1/TCOH
Definition: LUT.h:266
HOUGHptfLUT
This structure stores the patch-time-frequency look up table.
Definition: LUT.h:246
HOUGHptfLUT::offset
INT4 offset
Frequency bin corresponding to center of patch measured with respect to f0Bin (zero in modulated case...
Definition: LUT.h:255
HOUGHptfLUT::deltaF
REAL8 deltaF
Frequency resolution df=1/TCOH, where 1/TCOH is the coherent integration time used in teh demodulatio...
Definition: LUT.h:249
HOUGHptfLUT::maxNBins
UINT2 maxNBins
Maximum number of bins affecting the patch (for memory allocation purposes)
Definition: LUT.h:256
HOUGHptfLUT::nBin
INT4 nBin
Exact number of bins affecting the patch.
Definition: LUT.h:254
HOUGHptfLUT::border
HOUGHBorder * border
The annulus borders.
Definition: LUT.h:258
HOUGHptfLUT::iniBin
INT4 iniBin
First bin affecting the patch with respect to f0.
Definition: LUT.h:253
HOUGHptfLUT::bin
HOUGHBin2Border * bin
Bin to border correspondence.
Definition: LUT.h:259
HOUGHptfLUT::maxNBorders
UINT2 maxNBorders
Maximum number of borders affecting the patch (for memory allocation purposes)
Definition: LUT.h:257
HOUGHptfLUT::f0Bin
INT8 f0Bin
Frequency bin for which it has been constructed.
Definition: LUT.h:248
HOUGHptfLUT::nFreqValid
INT8 nFreqValid
Number of frequencies where the lut is valid.
Definition: LUT.h:252