Tesseract 3.02

- Main Page
- Related Pages
- Modules
- Namespaces
- Classes
- Files

- File List
- File Members

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
intfx.cpp
Go to the documentation of this file.
1 /******************************************************************************
2 ** Filename: intfx.c
3 ** Purpose: Integer character normalization & feature extraction
4 ** Author: Robert Moss
5 ** History: Tue May 21 15:51:57 MDT 1991, RWM, Created.
6 **
7 ** (c) Copyright Hewlett-Packard Company, 1988.
8 ** Licensed under the Apache License, Version 2.0 (the "License");
9 ** you may not use this file except in compliance with the License.
10 ** You may obtain a copy of the License at
11 ** http://www.apache.org/licenses/LICENSE-2.0
12 ** Unless required by applicable law or agreed to in writing, software
13 ** distributed under the License is distributed on an "AS IS" BASIS,
14 ** WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
15 ** See the License for the specific language governing permissions and
16 ** limitations under the License.
17 ******************************************************************************/
21 #include "intfx.h"
22 #include "intmatcher.h"
23 #include "const.h"
24 #include "helpers.h"
25 #include "ccutil.h"
26 #include "statistc.h"
27 #include "trainingsample.h"
28 #ifdef __UNIX__
29 #endif
30
31 using tesseract::TrainingSample;
32
36 int SaveFeature();
37 uinT8 BinaryAnglePlusPi(inT32 Y, inT32 X);
38 uinT8 MySqrt2();
39 void ClipRadius();
40
41 INT_VAR(classify_radius_gyr_min_man, 255,
42 "Minimum Radius of Gyration Mantissa 0-255: ");
43
44 INT_VAR(classify_radius_gyr_min_exp, 0,
45 "Minimum Radius of Gyration Exponent 0-255: ");
46
47 INT_VAR(classify_radius_gyr_max_man, 158,
48 "Maximum Radius of Gyration Mantissa 0-255: ");
49
50 INT_VAR(classify_radius_gyr_max_exp, 8,
51 "Maximum Radius of Gyration Exponent 0-255: ");
52
56 #define ATAN_TABLE_SIZE 64
57
58 // Look up table for arc tangent containing:
59 // atan(0.0) ... atan(ATAN_TABLE_SIZE - 1 / ATAN_TABLE_SIZE)
60 // The entries are in binary degrees where a full circle is 256 binary degrees.
61 static uinT8 AtanTable[ATAN_TABLE_SIZE];
62 // Look up table for cos and sin to turn the intfx feature angle to a vector.
63 // Also protected by atan_table_mutex.
64 static float cos_table[INT_CHAR_NORM_RANGE];
65 static float sin_table[INT_CHAR_NORM_RANGE];
66 // Guards write access to AtanTable so we dont create it more than once.
67 tesseract::CCUtilMutex atan_table_mutex;
68
69
73 /*---------------------------------------------------------------------------*/
74 void InitIntegerFX() {
75 static bool atan_table_init = false;
76 atan_table_mutex.Lock();
77 if (!atan_table_init) {
78 for (int i = 0; i < ATAN_TABLE_SIZE; i++) {
79 AtanTable[i] =
80 (uinT8) (atan ((i / (float) ATAN_TABLE_SIZE)) * 128.0 / PI + 0.5);
81 }
82 for (int i = 0; i < INT_CHAR_NORM_RANGE; ++i) {
83 cos_table[i] = cos(i * 2 * PI / INT_CHAR_NORM_RANGE + PI);
84 sin_table[i] = sin(i * 2 * PI / INT_CHAR_NORM_RANGE + PI);
85 }
86 atan_table_init = true;
87 }
88 atan_table_mutex.Unlock();
89 }
90
91 // Returns a vector representing the direction of a feature with the given
92 // theta direction in an INT_FEATURE_STRUCT.
93 FCOORD FeatureDirection(uinT8 theta) {
94 return FCOORD(cos_table[theta], sin_table[theta]);
95 }
96
97 TrainingSample* GetIntFeatures(tesseract::NormalizationMode mode,
98 TBLOB *blob, const DENORM& denorm) {
99 INT_FEATURE_ARRAY blfeatures;
100 INT_FEATURE_ARRAY cnfeatures;
101 INT_FX_RESULT_STRUCT fx_info;
102 ExtractIntFeat(blob, denorm, blfeatures, cnfeatures, &fx_info, NULL);
103 TrainingSample* sample = NULL;
104 if (mode == tesseract::NM_CHAR_ANISOTROPIC) {
105 int num_features = fx_info.NumCN;
106 if (num_features > 0) {
107 sample = TrainingSample::CopyFromFeatures(fx_info, cnfeatures,
108 num_features);
109 }
110 } else if (mode == tesseract::NM_BASELINE) {
111 int num_features = fx_info.NumBL;
112 if (num_features > 0) {
113 sample = TrainingSample::CopyFromFeatures(fx_info, blfeatures,
114 num_features);
115 }
116 } else {
117 ASSERT_HOST(!"Unsupported normalization mode!");
118 }
119 return sample;
120 }
121
122
123 /*--------------------------------------------------------------------------*/
124 // Extract a set of standard-sized features from Blobs and write them out in
125 // two formats: baseline normalized and character normalized.
126 //
127 // We presume the Blobs are already scaled so that x-height=128 units
128 //
129 // Standard Features:
130 // We take all outline segments longer than 7 units and chop them into
131 // standard-sized segments of approximately 13 = (64 / 5) units.
132 // When writing these features out, we output their center and angle as
133 // measured counterclockwise from the vector <-1, 0>
134 //
135 // Baseline Normalized Output:
136 // We center the grapheme by aligning the x-coordinate of its centroid with
137 // x=0 and subtracting 128 from the y-coordinate.
138 //
139 // Character Normalized Output:
140 // We align the grapheme's centroid at the origin and scale it asymmetrically
141 // in x and y so that the result is vaguely square.
142 //
143 int ExtractIntFeat(TBLOB *Blob,
144 const DENORM& denorm,
145 INT_FEATURE_ARRAY BLFeat,
146 INT_FEATURE_ARRAY CNFeat,
147 INT_FX_RESULT_STRUCT* Results,
148 inT32 *FeatureOutlineArray) {
149
150 TESSLINE *OutLine;
151 EDGEPT *Loop, *LoopStart, *Segment;
152 inT16 LastX, LastY, Xmean, Ymean;
153 inT32 NormX, NormY, DeltaX, DeltaY;
154 inT32 Xsum, Ysum;
155 uinT32 Ix, Iy, LengthSum;
156 uinT16 n;
157 // n - the number of features to extract from a given outline segment.
158 // We extract features from every outline segment longer than ~6 units.
159 // We chop these long segments into standard-sized features approximately
160 // 13 (= 64 / 5) units in length.
161 uinT8 Theta;
162 uinT16 NumBLFeatures, NumCNFeatures;
163 uinT8 RxInv, RyInv; /* x.xxxxxxx * 2^Exp */
164 uinT8 RxExp, RyExp;
165 /* sxxxxxxxxxxxxxxxxxxxxxxx.xxxxxxxx */
166 register inT32 pfX, pfY, dX, dY;
167 uinT16 Length;
168 register int i;
169
170 Results->Length = 0;
171 Results->Xmean = 0;
172 Results->Ymean = 0;
173 Results->Rx = 0;
174 Results->Ry = 0;
175 Results->NumBL = 0;
176 Results->NumCN = 0;
177 Results->YBottom = MAX_UINT8;
178 Results->YTop = 0;
179
180 // Calculate the centroid (Xmean, Ymean) for the blob.
181 // We use centroid (instead of center of bounding box or center of smallest
182 // enclosing circle) so the algorithm will not be too greatly influenced by
183 // small amounts of information at the edge of a character's bounding box.
184 NumBLFeatures = 0;
185 NumCNFeatures = 0;
186 OutLine = Blob->outlines;
187 Xsum = 0;
188 Ysum = 0;
189 LengthSum = 0;
190 while (OutLine != NULL) {
191 LoopStart = OutLine->loop;
192 Loop = LoopStart;
193 LastX = Loop->pos.x;
194 LastY = Loop->pos.y;
195 /* Check for bad loops */
196 if ((Loop == NULL) || (Loop->next == NULL) || (Loop->next == LoopStart))
197 return FALSE;
198 do {
199 Segment = Loop;
200 Loop = Loop->next;
201 NormX = Loop->pos.x;
202 NormY = Loop->pos.y;
203
204 n = 1;
205 if (!Segment->IsHidden()) {
206 DeltaX = NormX - LastX;
207 DeltaY = NormY - LastY;
208 Length = MySqrt(DeltaX, DeltaY);
209 n = ((Length << 2) + Length + 32) >> 6;
210 if (n != 0) {
211 Xsum += ((LastX << 1) + DeltaX) * (int) Length;
212 Ysum += ((LastY << 1) + DeltaY) * (int) Length;
213 LengthSum += Length;
214 }
215 }
216 if (n != 0) { /* Throw away a point that is too close */
217 LastX = NormX;
218 LastY = NormY;
219 }
220 }
221 while (Loop != LoopStart);
222 OutLine = OutLine->next;
223 }
224 if (LengthSum == 0)
225 return FALSE;
226 Xmean = (Xsum / (inT32) LengthSum) >> 1;
227 Ymean = (Ysum / (inT32) LengthSum) >> 1;
228
229 Results->Length = LengthSum;
230 Results->Xmean = Xmean;
231 Results->Ymean = Ymean;
232
233 // Extract Baseline normalized features,
234 // and find 2nd moments (Ix, Iy) & radius of gyration (Rx, Ry).
235 //
236 // Ix = Sum y^2 dA, where:
237 // Ix: the second moment of area about the axis x
238 // dA = 1 for our standard-sized piece of outline
239 // y: the perependicular distance to the x axis
240 // Rx = sqrt(Ix / A)
241 // Note: 1 <= Rx <= height of blob / 2
242 // Ry = sqrt(Iy / A)
243 // Note: 1 <= Ry <= width of blob / 2
244 Ix = 0;
245 Iy = 0;
246 NumBLFeatures = 0;
247 OutLine = Blob->outlines;
248 int min_x = 0;
249 int max_x = 0;
250 while (OutLine != NULL) {
251 LoopStart = OutLine->loop;
252 Loop = LoopStart;
253 LastX = Loop->pos.x - Xmean;
254 LastY = Loop->pos.y;
255 /* Check for bad loops */
256 if ((Loop == NULL) || (Loop->next == NULL) || (Loop->next == LoopStart))
257 return FALSE;
258 do {
259 Segment = Loop;
260 Loop = Loop->next;
261 NormX = Loop->pos.x - Xmean;
262 NormY = Loop->pos.y;
263 if (NormY < Results->YBottom)
264 Results->YBottom = ClipToRange(NormY, 0, MAX_UINT8);
265 if (NormY > Results->YTop)
266 Results->YTop = ClipToRange(NormY, 0, MAX_UINT8);
267 UpdateRange(NormX, &min_x, &max_x);
268
269 n = 1;
270 if (!Segment->IsHidden()) {
271 DeltaX = NormX - LastX;
272 DeltaY = NormY - LastY;
273 Length = MySqrt(DeltaX, DeltaY);
274 n = ((Length << 2) + Length + 32) >> 6;
275 if (n != 0) {
276 Theta = BinaryAnglePlusPi(DeltaY, DeltaX);
277 dX = (DeltaX << 8) / n;
278 dY = (DeltaY << 8) / n;
279 pfX = (LastX << 8) + (dX >> 1);
280 pfY = (LastY << 8) + (dY >> 1);
281 Ix += ((pfY >> 8) - Ymean) * ((pfY >> 8) - Ymean);
282 // TODO(eger): Hmmm... Xmean is not necessarily 0.
283 // Figure out if we should center against Xmean for these
284 // features, and if so fix Iy & SaveFeature().
285 Iy += (pfX >> 8) * (pfX >> 8);
286 if (SaveFeature(BLFeat,
287 NumBLFeatures,
288 (inT16) (pfX >> 8),
289 (inT16) ((pfY >> 8) - 128),
290 Theta) == FALSE)
291 return FALSE;
292 NumBLFeatures++;
293 for (i = 1; i < n; i++) {
294 pfX += dX;
295 pfY += dY;
296 Ix += ((pfY >> 8) - Ymean) * ((pfY >> 8) - Ymean);
297 Iy += (pfX >> 8) * (pfX >> 8);
298 if (SaveFeature(BLFeat,
299 NumBLFeatures,
300 (inT16) (pfX >> 8),
301 (inT16) ((pfY >> 8) - 128),
302 Theta) == FALSE)
303 return FALSE;
304 NumBLFeatures++;
305 }
306 }
307 }
308 if (n != 0) { /* Throw away a point that is too close */
309 LastX = NormX;
310 LastY = NormY;
311 }
312 }
313 while (Loop != LoopStart);
314 OutLine = OutLine->next;
315 }
316 Results->Width = max_x - min_x;
317 if (Ix == 0)
318 Ix = 1;
319 if (Iy == 0)
320 Iy = 1;
321 RxInv = MySqrt2 (NumBLFeatures, Ix, &RxExp);
322 RyInv = MySqrt2 (NumBLFeatures, Iy, &RyExp);
323 ClipRadius(&RxInv, &RxExp, &RyInv, &RyExp);
324
325 Results->Rx = (inT16) (51.2 / (double) RxInv * pow (2.0, (double) RxExp));
326 Results->Ry = (inT16) (51.2 / (double) RyInv * pow (2.0, (double) RyExp));
327 if (Results->Ry == 0) {
328 /*
329 This would result in features having 'nan' values.
330 Since the expression is always > 0, assign a value of 1.
331 */
332 Results->Ry = 1;
333 }
334 Results->NumBL = NumBLFeatures;
335
336 // Extract character normalized features
337 //
338 // Rescale the co-ordinates to "equalize" distribution in X and Y, making
339 // all of the following unichars be sized to look similar: , ' 1 i
340 //
341 // We calculate co-ordinates relative to the centroid, and then scale them
342 // as follows (accomplishing a scale of up to 102.4 / dimension):
343 // y *= 51.2 / Rx [ y scaled by 0.0 ... 102.4 / height of glyph ]
344 // x *= 51.2 / Ry [ x scaled by 0.0 ... 102.4 / width of glyph ]
345 // Although tempting to think so, this does not guarantee that our range
346 // is within [-102.4...102.4] x [-102.4...102.4] because (Xmean, Ymean)
347 // is the centroid, not the center of the bounding box. Instead, we can
348 // only bound the result to [-204 ... 204] x [-204 ... 204]
349 //
350 NumCNFeatures = 0;
351 OutLine = Blob->outlines;
352 int OutLineIndex = -1;
353 while (OutLine != NULL) {
354 LoopStart = OutLine->loop;
355 Loop = LoopStart;
356 LastX = (Loop->pos.x - Xmean) * RyInv;
357 LastY = (Loop->pos.y - Ymean) * RxInv;
358 LastX >>= (inT8) RyExp;
359 LastY >>= (inT8) RxExp;
360 OutLineIndex++;
361
362 /* Check for bad loops */
363 if ((Loop == NULL) || (Loop->next == NULL) || (Loop->next == LoopStart))
364 return FALSE;
365 do {
366 Segment = Loop;
367 Loop = Loop->next;
368 NormX = (Loop->pos.x - Xmean) * RyInv;
369 NormY = (Loop->pos.y - Ymean) * RxInv;
370 NormX >>= (inT8) RyExp;
371 NormY >>= (inT8) RxExp;
372
373 n = 1;
374 if (!Segment->IsHidden()) {
375 DeltaX = NormX - LastX;
376 DeltaY = NormY - LastY;
377 Length = MySqrt(DeltaX, DeltaY);
378 n = ((Length << 2) + Length + 32) >> 6;
379 if (n != 0) {
380 Theta = BinaryAnglePlusPi(DeltaY, DeltaX);
381 dX = (DeltaX << 8) / n;
382 dY = (DeltaY << 8) / n;
383 pfX = (LastX << 8) + (dX >> 1);
384 pfY = (LastY << 8) + (dY >> 1);
385 if (SaveFeature(CNFeat,
386 NumCNFeatures,
387 (inT16) (pfX >> 8),
388 (inT16) (pfY >> 8),
389 Theta) == FALSE)
390 return FALSE;
391 if (FeatureOutlineArray) {
392 FeatureOutlineArray[NumCNFeatures] = OutLineIndex;
393 }
394 NumCNFeatures++;
395 for (i = 1; i < n; i++) {
396 pfX += dX;
397 pfY += dY;
398 if (SaveFeature(CNFeat,
399 NumCNFeatures,
400 (inT16) (pfX >> 8),
401 (inT16) (pfY >> 8),
402 Theta) == FALSE)
403 return FALSE;
404 if (FeatureOutlineArray) {
405 FeatureOutlineArray[NumCNFeatures] = OutLineIndex;
406 }
407 NumCNFeatures++;
408 }
409 }
410 }
411 if (n != 0) { /* Throw away a point that is too close */
412 LastX = NormX;
413 LastY = NormY;
414 }
415 }
416 while (Loop != LoopStart);
417 OutLine = OutLine->next;
418 }
419
420 Results->NumCN = NumCNFeatures;
421 return TRUE;
422 }
423
424
425 /*--------------------------------------------------------------------------*/
426 // Return the "binary angle" [0..255]
427 // made by vector <X, Y> as measured counterclockwise from <-1, 0>
428 // The order of the arguments follows the convention of atan2(3)
429 uinT8 BinaryAnglePlusPi(inT32 Y, inT32 X) {
430 inT16 Angle, Atan;
431 uinT16 Ratio;
432 uinT32 AbsX, AbsY;
433
434 assert ((X != 0) || (Y != 0));
435 if (X < 0)
436 AbsX = -X;
437 else
438 AbsX = X;
439 if (Y < 0)
440 AbsY = -Y;
441 else
442 AbsY = Y;
443 if (AbsX > AbsY)
444 Ratio = AbsY * ATAN_TABLE_SIZE / AbsX;
445 else
446 Ratio = AbsX * ATAN_TABLE_SIZE / AbsY;
447 if (Ratio >= ATAN_TABLE_SIZE)
448 Ratio = ATAN_TABLE_SIZE - 1;
449 Atan = AtanTable[Ratio];
450 if (X >= 0)
451 if (Y >= 0)
452 if (AbsX > AbsY)
453 Angle = Atan;
454 else
455 Angle = 64 - Atan;
456 else if (AbsX > AbsY)
457 Angle = 256 - Atan;
458 else
459 Angle = 192 + Atan;
460 else if (Y >= 0)
461 if (AbsX > AbsY)
462 Angle = 128 - Atan;
463 else
464 Angle = 64 + Atan;
465 else if (AbsX > AbsY)
466 Angle = 128 + Atan;
467 else
468 Angle = 192 - Atan;
469
470 /* reverse angles to match old feature extractor: Angle += PI */
471 Angle += 128;
472 Angle &= 255;
473 return (uinT8) Angle;
474 }
475
476
477 /*--------------------------------------------------------------------------*/
478 int SaveFeature(INT_FEATURE_ARRAY FeatureArray,
479 uinT16 FeatureNum,
480 inT16 X,
481 inT16 Y,
482 uinT8 Theta) {
483 INT_FEATURE Feature;
484
485 if (FeatureNum >= MAX_NUM_INT_FEATURES)
486 return FALSE;
487
488 Feature = &(FeatureArray[FeatureNum]);
489
490 X = X + 128;
491 Y = Y + 128;
492
493 Feature->X = ClipToRange<inT16>(X, 0, 255);
494 Feature->Y = ClipToRange<inT16>(Y, 0, 255);
495 Feature->Theta = Theta;
496 Feature->CP_misses = 0;
497
498 return TRUE;
499 }
500
501
502 /*---------------------------------------------------------------------------*/
503 // Return floor(sqrt(min(emm, x)^2 + min(emm, y)^2))
504 // where emm = EvidenceMultMask.
505 uinT16 MySqrt(inT32 X, inT32 Y) {
506 register uinT16 SqRoot;
507 register uinT32 Square;
508 register uinT16 BitLocation;
509 register uinT32 Sum;
510 const uinT32 EvidenceMultMask =
511 ((1 << IntegerMatcher::kIntEvidenceTruncBits) - 1);
512
513 if (X < 0)
514 X = -X;
515 if (Y < 0)
516 Y = -Y;
517
518 if (X > EvidenceMultMask)
519 X = EvidenceMultMask;
520 if (Y > EvidenceMultMask)
521 Y = EvidenceMultMask;
522
523 Sum = X * X + Y * Y;
524
525 BitLocation = (EvidenceMultMask + 1) << 1;
526 SqRoot = 0;
527 do {
528 Square = (SqRoot | BitLocation) * (SqRoot | BitLocation);
529 if (Square <= Sum)
530 SqRoot |= BitLocation;
531 BitLocation >>= 1;
532 }
533 while (BitLocation);
534
535 return SqRoot;
536 }
537
538
539 /*--------------------------------------------------------------------------*/
540 // Return two integers which can be used to express the sqrt(I/N):
541 // sqrt(I/N) = 51.2 * 2^(*Exp) / retval
542 uinT8 MySqrt2(uinT16 N, uinT32 I, uinT8 *Exp) {
543 register inT8 k;
544 register uinT32 N2;
545 register uinT8 SqRoot;
546 register uinT16 Square;
547 register uinT8 BitLocation;
548 register uinT16 Ratio;
549
550 N2 = N * 41943;
551
552 k = 9;
553 while ((N2 & 0xc0000000) == 0) {
554 N2 <<= 2;
555 k += 1;
556 }
557
558 while ((I & 0xc0000000) == 0) {
559 I <<= 2;
560 k -= 1;
561 }
562
563 if (((N2 & 0x80000000) == 0) && ((I & 0x80000000) == 0)) {
564 N2 <<= 1;
565 I <<= 1;
566 }
567
568 N2 &= 0xffff0000;
569 I >>= 14;
570 Ratio = N2 / I;
571
572 BitLocation = 128;
573 SqRoot = 0;
574 do {
575 Square = (SqRoot | BitLocation) * (SqRoot | BitLocation);
576 if (Square <= Ratio)
577 SqRoot |= BitLocation;
578 BitLocation >>= 1;
579 }
580 while (BitLocation);
581
582 if (k < 0) {
583 *Exp = 0;
584 return 255;
585 }
586 else {
587 *Exp = k;
588 return SqRoot;
589 }
590 }
591
592
593 /*-------------------------------------------------------------------------*/
594 void ClipRadius(uinT8 *RxInv, uinT8 *RxExp, uinT8 *RyInv, uinT8 *RyExp) {
595 register uinT8 AM, BM, AE, BE;
596 register uinT8 BitN, LastCarry;
597 int RxInvLarge, RyInvSmall;
598
599 AM = classify_radius_gyr_min_man;
600 AE = classify_radius_gyr_min_exp;
601 BM = *RxInv;
602 BE = *RxExp;
603 LastCarry = 1;
604 while ((AM != 0) || (BM != 0)) {
605 if (AE > BE) {
606 BitN = LastCarry + (AM & 1) + 1;
607 AM >>= 1;
608 AE--;
609 }
610 else if (AE < BE) {
611 BitN = LastCarry + (!(BM & 1));
612 BM >>= 1;
613 BE--;
614 }
615 else { /* AE == BE */
616 BitN = LastCarry + (AM & 1) + (!(BM & 1));
617 AM >>= 1;
618 BM >>= 1;
619 AE--;
620 BE--;
621 }
622 LastCarry = (BitN & 2) > 1;
623 BitN = BitN & 1;
624 }
625 BitN = LastCarry + 1;
626 LastCarry = (BitN & 2) > 1;
627 BitN = BitN & 1;
628
629 if (BitN == 1) {
630 *RxInv = classify_radius_gyr_min_man;
631 *RxExp = classify_radius_gyr_min_exp;
632 }
633
634 AM = classify_radius_gyr_min_man;
635 AE = classify_radius_gyr_min_exp;
636 BM = *RyInv;
637 BE = *RyExp;
638 LastCarry = 1;
639 while ((AM != 0) || (BM != 0)) {
640 if (AE > BE) {
641 BitN = LastCarry + (AM & 1) + 1;
642 AM >>= 1;
643 AE--;
644 }
645 else if (AE < BE) {
646 BitN = LastCarry + (!(BM & 1));
647 BM >>= 1;
648 BE--;
649 }
650 else { /* AE == BE */
651 BitN = LastCarry + (AM & 1) + (!(BM & 1));
652 AM >>= 1;
653 BM >>= 1;
654 AE--;
655 BE--;
656 }
657 LastCarry = (BitN & 2) > 1;
658 BitN = BitN & 1;
659 }
660 BitN = LastCarry + 1;
661 LastCarry = (BitN & 2) > 1;
662 BitN = BitN & 1;
663
664 if (BitN == 1) {
665 *RyInv = classify_radius_gyr_min_man;
666 *RyExp = classify_radius_gyr_min_exp;
667 }
668
669 AM = classify_radius_gyr_max_man;
670 AE = classify_radius_gyr_max_exp;
671 BM = *RxInv;
672 BE = *RxExp;
673 LastCarry = 1;
674 while ((AM != 0) || (BM != 0)) {
675 if (AE > BE) {
676 BitN = LastCarry + (AM & 1) + 1;
677 AM >>= 1;
678 AE--;
679 }
680 else if (AE < BE) {
681 BitN = LastCarry + (!(BM & 1));
682 BM >>= 1;
683 BE--;
684 }
685 else { /* AE == BE */
686 BitN = LastCarry + (AM & 1) + (!(BM & 1));
687 AM >>= 1;
688 BM >>= 1;
689 AE--;
690 BE--;
691 }
692 LastCarry = (BitN & 2) > 1;
693 BitN = BitN & 1;
694 }
695 BitN = LastCarry + 1;
696 LastCarry = (BitN & 2) > 1;
697 BitN = BitN & 1;
698
699 if (BitN == 1)
700 RxInvLarge = 1;
701 else
702 RxInvLarge = 0;
703
704 AM = *RyInv;
705 AE = *RyExp;
706 BM = classify_radius_gyr_max_man;
707 BE = classify_radius_gyr_max_exp;
708 LastCarry = 1;
709 while ((AM != 0) || (BM != 0)) {
710 if (AE > BE) {
711 BitN = LastCarry + (AM & 1) + 1;
712 AM >>= 1;
713 AE--;
714 }
715 else if (AE < BE) {
716 BitN = LastCarry + (!(BM & 1));
717 BM >>= 1;
718 BE--;
719 }
720 else { /* AE == BE */
721 BitN = LastCarry + (AM & 1) + (!(BM & 1));
722 AM >>= 1;
723 BM >>= 1;
724 AE--;
725 BE--;
726 }
727 LastCarry = (BitN & 2) > 1;
728 BitN = BitN & 1;
729 }
730 BitN = LastCarry + 1;
731 LastCarry = (BitN & 2) > 1;
732 BitN = BitN & 1;
733
734 if (BitN == 1)
735 RyInvSmall = 1;
736 else
737 RyInvSmall = 0;
738
739 if (RxInvLarge && RyInvSmall) {
740 *RyInv = classify_radius_gyr_max_man;
741 *RyExp = classify_radius_gyr_max_exp;
742 }
743
744 }