Loading [MathJax]/extensions/TeX/AMSsymbols.js
LALPulsar 7.1.1.1-5e288d3
•All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Modules Pages
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
98static void PLUTInitialize( HOUGHptfLUT * );
99static void FillPLUT( HOUGHParamPLUT *, HOUGHptfLUT *, HOUGHPatchGrid * );
100
101static void CheckLeftCircle( REAL8, REAL8, REAL8, INT4 *, INT4 *, INT4 *,
102 HOUGHPatchGrid * );
103static void CheckRightCircle( REAL8, REAL8, REAL8, INT4 *, INT4 *, INT4 *,
104 HOUGHPatchGrid * );
105static void DrawRightCircle( REAL8, REAL8, REAL8, INT4, INT4, COORType *,
106 HOUGHPatchGrid * );
107static void DrawLeftCircle( REAL8, REAL8, REAL8, INT4, INT4, COORType *,
108 HOUGHPatchGrid * );
109static void CheckLineCase( REAL8, REAL8, REAL8, REAL8 *, INT4 * );
110/* static void FindExactLine(REAL8,REAL8,REAL8 *,REAL8 *); */
111static void FindLine( REAL8, REAL8, REAL8, REAL8 *, REAL8 * );
112static void CheckLineIntersection( REAL8, REAL8, REAL8, INT4 *, INT4 *, INT4 *,
113 HOUGHPatchGrid * );
114static void DrawLine( REAL8, REAL8, REAL8, INT4, INT4, COORType *, HOUGHPatchGrid * );
115static void Fill1Column( INT4, UINT4 *, HOUGHptfLUT *, HOUGHPatchGrid * );
116static void FillCaseN1( INT4, INT4, INT4, HOUGHptfLUT *, HOUGHPatchGrid * );
117static void FillCaseN2( INT4, INT4, HOUGHptfLUT *, HOUGHPatchGrid * );
118static void FillCaseN3( INT4, INT4, INT4, INT4 *, HOUGHptfLUT *, HOUGHPatchGrid * );
119static void FillCaseN4( INT4, INT4, HOUGHptfLUT *, HOUGHPatchGrid * );
120static void FillCaseN5( INT4, INT4, INT4, HOUGHptfLUT * );
121static void FillCaseN6( INT4, INT4, INT4, HOUGHptfLUT *, HOUGHPatchGrid * );
122static void FillCaseN7( INT4, INT4, HOUGHptfLUT *, HOUGHPatchGrid * );
123static void FillCaseN8( INT4, INT4, HOUGHptfLUT *, HOUGHPatchGrid * );
124static void Fill1ColumnAnor( INT4, HOUGHptfLUT *, HOUGHPatchGrid * );
125static void FillCaseA1( INT4, INT4, INT4, HOUGHptfLUT * );
126static void FillCaseA2( INT4, INT4, INT4, HOUGHptfLUT * );
127static void FillCaseA3( INT4, INT4, INT4, HOUGHptfLUT *, HOUGHPatchGrid * );
128
129static void InitialCircleCase( UINT4 *, REAL8, REAL8, REAL8, REAL8 *, INT4 *, INT4 *,
130 HOUGHptfLUT *, HOUGHPatchGrid * );
131static void SecondCircleCase( INT4, UINT4 *, REAL8, REAL8, REAL8, INT4, REAL8 *,
132 INT4 *, INT4 *, INT4 *, HOUGHptfLUT *, HOUGHPatchGrid * );
133static void FollowCircleCase( INT4, UINT4 *, REAL8, REAL8, REAL8, REAL8, REAL8, INT4 *,
134 INT4 *, INT4 *, HOUGHptfLUT *, HOUGHPatchGrid * );
135static void InitialLineCase( UINT4 *, REAL8, REAL8, REAL8, INT4 *, HOUGHptfLUT *,
136 HOUGHPatchGrid * );
137static void SecondLineCase( INT4, UINT4 *, REAL8, REAL8, REAL8, INT4 *, HOUGHptfLUT *,
138 HOUGHPatchGrid * );
139static 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/* >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< */
149void 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
244static 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
332static 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
626static 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
759static 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
899static 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
1105static 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 */
1164static 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
1189static 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
1213static 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/*****************************************************************/
1291static 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
1374static 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
1612static 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
1868static 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
2128static 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
2340static 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
2546static 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
2586static 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
2639static 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/***************************************************************/
2723static 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/***************************************************************/
2768static 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/***************************************************************/
2794static 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/***************************************************************/
2842static 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/***************************************************************/
2873static 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/***************************************************************/
2908static 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/***************************************************************/
2951static 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/***************************************************************/
2977static 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
3023static 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/***************************************************************/
3107static 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/***************************************************************/
3138static 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/***************************************************************/
3169static 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:2446
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