ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
ptwXY_methods.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file ptwXY_methods.cc
1 /*
2 # <<BEGIN-copyright>>
3 # <<END-copyright>>
4 */
5 
6 #include <cmath>
7 #include <float.h>
8 
9 #include "ptwXY.h"
10 
11 #if defined __cplusplus
12 #include <cmath>
13 #include "G4Exp.hh"
14 #include "G4Log.hh"
15 namespace GIDI {
16 using namespace GIDI;
17 #endif
18 
19 static nfu_status ptwXY_clip2( ptwXYPoints *ptwXY1, double y, double x1, double y1, double x2, double y2 );
20 static double ptwXY_thicken_linear_dx( int sectionSubdivideMax, double dxMax, double x1, double x2 );
21 static nfu_status ptwXY_thin2( ptwXYPoints *thinned, char *thin, double accuracy, int64_t i1, int64_t i2 );
22 /*
23 ************************************************************
24 */
25 nfu_status ptwXY_clip( ptwXYPoints *ptwXY1, double yMin, double yMax ) {
26 /*
27  This function acts oddly for xy = [ [ 1, 0 ], [ 3, -2 ], [ 4, 1 ] ] and yMin = 0.2, why???????
28  This function probably only works for linear, linear interpolation (mainly because of ptwXY_clip2).
29 */
30  int64_t i, j, n;
31  double x2, y2;
32  nfu_status status;
33  ptwXYPoints *clipped;
34  ptwXYPoint *points;
35 
36  if( ( status = ptwXY_simpleCoalescePoints( ptwXY1 ) ) != nfu_Okay ) return( status );
38  n = ptwXY1->length;
39  if( n > 0 ) {
40  i = 0;
41  if( ptwXY_getYMax( ptwXY1 ) < yMin ) i = 1;
42  if( ptwXY_getYMin( ptwXY1 ) > yMax ) i = 1;
43  if( i == 1 ) return( ptwXY_clear( ptwXY1 ) );
44  }
45  if( n == 1 ) {
46  y2 = ptwXY1->points[0].y;
47  if( y2 < yMin ) {
48  ptwXY1->points[0].y = yMin; }
49  else if( y2 > yMax ) {
50  ptwXY1->points[0].y = yMax;
51  } }
52  else if( n > 1 ) {
53  if( ( clipped = ptwXY_new( ptwXY1->interpolation, &(ptwXY1->interpolationOtherInfo),
54  ptwXY1->biSectionMax, ptwXY1->accuracy, n, 10, &status, ptwXY1->userFlag ) ) == NULL )
55  return( ptwXY1->status = status );
56  for( i = 0; i < n; i++ ) {
57  x2 = ptwXY1->points[i].x;
58  y2 = ptwXY1->points[i].y;
59  if( y2 < yMin ) {
60  if( i > 0 ) {
61  points = ptwXY_getPointAtIndex_Unsafely( clipped, clipped->length - 1 );
62  if( points->y > yMin ) {
63  if( ( status = ptwXY_clip2( clipped, yMin, points->x, points->y, x2, y2 ) ) != nfu_Okay ) goto Err;
64  }
65  }
66  if( ( status = ptwXY_setValueAtX( clipped, x2, yMin ) ) != nfu_Okay ) goto Err;
67  j = i;
68  for( i++; i < n; i++ ) if( !( ptwXY1->points[i].y < yMin ) ) break;
69  if( i < n ) {
70  x2 = ptwXY1->points[i].x;
71  y2 = ptwXY1->points[i].y;
72  if( ( status = ptwXY_clip2( clipped, yMin, ptwXY1->points[i-1].x, ptwXY1->points[i-1].y, x2, y2 ) ) != nfu_Okay ) goto Err;
73  if( y2 > yMax ) {
74  if( ( status = ptwXY_clip2( clipped, yMax, ptwXY1->points[i-1].x, ptwXY1->points[i-1].y, x2, y2 ) ) != nfu_Okay ) goto Err;
75  } }
76  else if( j != n - 1 ) {
77  if( ( status = ptwXY_setValueAtX( clipped, ptwXY1->points[n - 1].x, yMin ) ) != nfu_Okay ) goto Err;
78  }
79  i--; }
80  else if( y2 > yMax ) {
81  if( i > 0 ) {
82  points = ptwXY_getPointAtIndex_Unsafely( clipped, clipped->length - 1 );
83  if( points->y < yMax ) {
84  if( ( status = ptwXY_clip2( clipped, yMax, points->x, points->y, x2, y2 ) ) != nfu_Okay ) goto Err;
85  }
86  }
87  if( ( status = ptwXY_setValueAtX( clipped, x2, yMax ) ) != nfu_Okay ) goto Err;
88  j = i;
89  for( i++; i < n; i++ ) if( !( ptwXY1->points[i].y > yMax ) ) break;
90  if( i < n ) {
91  x2 = ptwXY1->points[i].x;
92  y2 = ptwXY1->points[i].y;
93  if( ( status = ptwXY_clip2( clipped, yMax, ptwXY1->points[i-1].x, ptwXY1->points[i-1].y, x2, y2 ) ) != nfu_Okay ) goto Err;
94  if( y2 < yMin ) {
95  if( ( status = ptwXY_clip2( clipped, yMin, ptwXY1->points[i-1].x, ptwXY1->points[i-1].y, x2, y2 ) ) != nfu_Okay ) goto Err;
96  } }
97  else if( j != n - 1 ) {
98  if( ( status = ptwXY_setValueAtX( clipped, ptwXY1->points[n - 1].x, yMax ) ) != nfu_Okay ) goto Err;
99  }
100  i--; }
101  else {
102  if( ( status = ptwXY_setValueAtX( clipped, x2, y2 ) ) != nfu_Okay ) goto Err;
103  }
104  }
105  if( ( status = ptwXY_simpleCoalescePoints( clipped ) ) != nfu_Okay ) goto Err;
106  ptwXY1->length = clipped->length; /* The squeamish may want to skip the next few lines. */
107  clipped->length = n;
108  n = ptwXY1->allocatedSize;
109  ptwXY1->allocatedSize = clipped->allocatedSize;
110  clipped->allocatedSize = n;
111  points = clipped->points;
112  clipped->points = ptwXY1->points;
113  ptwXY1->points = points;
114  ptwXY_free( clipped );
115  }
116 
117  return( ptwXY1->status );
118 
119 Err:
120  ptwXY_free( clipped );
121  return( ptwXY1->status = status );
122 }
123 /*
124 ************************************************************
125 */
126 static nfu_status ptwXY_clip2( ptwXYPoints *clipped, double y, double x1, double y1, double x2, double y2 ) {
127 
128  double x;
129  nfu_status status = nfu_Okay;
130 
131  x = ( y - y1 ) * ( x2 - x1 ) / ( y2 - y1 ) + x1;
132  if( x <= x1 ) {
133  x = x1; }
134  else if( x >= x2 ) {
135  x = x1; }
136  else {
137  status = ptwXY_setValueAtX( clipped, x, y );
138  }
139  return( status );
140 }
141 /*
142 ************************************************************
143 */
144 nfu_status ptwXY_thicken( ptwXYPoints *ptwXY1, int sectionSubdivideMax, double dxMax, double fxMax ) {
145 
146  double x1, x2 = 0., y1, y2 = 0., fx = 1.1, x, dx, dxp, lfx, y; /* fx initialized so compilers want complain. */
147  int64_t i, notFirstPass = 0;
148  int nfx, nDone, doLinear;
149  nfu_status status;
150 
152  if( ( sectionSubdivideMax < 1 ) || ( dxMax < 0. ) || ( fxMax < 1. ) ) return( nfu_badInput );
153  if( sectionSubdivideMax > ptwXY_sectionSubdivideMax ) sectionSubdivideMax = ptwXY_sectionSubdivideMax;
154  if( ( status = ptwXY_simpleCoalescePoints( ptwXY1 ) ) != nfu_Okay ) return( status );
155  for( i = ptwXY1->length - 1; i >= 0; i-- ) {
156  x1 = ptwXY1->points[i].x;
157  y1 = ptwXY1->points[i].y;
158  if( notFirstPass ) {
159  dx = ptwXY_thicken_linear_dx( sectionSubdivideMax, dxMax, x1, x2 );
160 
161  if( x1 == 0. ) {
162  doLinear = 1; }
163  else {
164  fx = x2 / x1;
165  if( fx > 0. ) {
166  lfx = G4Log( fx );
167  if( fxMax == 1. ) {
168  nfx = sectionSubdivideMax; }
169  else {
170  nfx = ( (int) ( lfx / G4Log( fxMax ) ) ) + 1;
171  if( nfx > sectionSubdivideMax ) nfx = sectionSubdivideMax;
172  }
173  if( nfx > 0 ) fx = G4Exp( lfx / nfx );
174  doLinear = 0;
175  if( dx < ( fx - 1 ) * x1 ) doLinear = 1; }
176  else {
177  doLinear = 1;
178  }
179  }
180  x = x1;
181  dxp = dx;
182  nDone = 0;
183  while( 1 ) {
184  if( doLinear ) {
185  x += dx; }
186  else {
187  dx = ptwXY_thicken_linear_dx( sectionSubdivideMax - nDone, dxMax, x, x2 );
188  if( dx <= ( fx - 1 ) * x ) {
189  dxp = dx;
190  doLinear = 1;
191  continue;
192  }
193  dxp = ( fx - 1. ) * x;
194  x *= fx;
195  }
196  if( ( x2 - x ) < 0.05 * std::fabs( dxp ) ) break;
197  if( ( status = ptwXY_interpolatePoint( ptwXY1->interpolation, x, &y, x1, y1, x2, y2 ) ) != nfu_Okay ) return( status );
198  if( ( status = ptwXY_setValueAtX( ptwXY1, x, y ) ) != nfu_Okay ) return( status );
199  nDone++;
200  } // Loop checking, 11.06.2015, T. Koi
201  }
202  notFirstPass = 1;
203  x2 = x1;
204  y2 = y1;
205  }
206  return( status );
207 }
208 /*
209 ************************************************************
210 */
211 ptwXYPoints *ptwXY_thin( ptwXYPoints *ptwXY1, double accuracy, nfu_status *status ) {
212 
213  int64_t i, j, length = ptwXY1->length;
214  ptwXYPoints *thinned = NULL;
215  double y1, y2, y3;
216  char *thin = NULL;
217 
218  if( length < 3 ) return( ptwXY_clone( ptwXY1, status ) ); /* Logic below requires at least 2 points. */
219  if( ( *status = ptwXY_simpleCoalescePoints( ptwXY1 ) ) != nfu_Okay ) return( NULL );
220  *status = nfu_otherInterpolation;
221  if( ptwXY1->interpolation == ptwXY_interpolationOther ) return( NULL );
222  if( accuracy < ptwXY1->accuracy ) accuracy = ptwXY1->accuracy;
223  if( ( thinned = ptwXY_new( ptwXY1->interpolation, &(ptwXY1->interpolationOtherInfo),
224  ptwXY1->biSectionMax, accuracy, length, ptwXY1->overflowLength, status, ptwXY1->userFlag ) ) == NULL ) return( NULL );
225 
226  thinned->points[0] = ptwXY1->points[0]; /* This sections removes middle point if surrounding points have the same y-value. */
227  y1 = ptwXY1->points[0].y;
228  y2 = ptwXY1->points[1].y;
229  for( i = 2, j = 1; i < length; i++ ) {
230  y3 = ptwXY1->points[i].y;
231  if( ( y1 != y2 ) || ( y2 != y3 ) ) {
232  thinned->points[j++] = ptwXY1->points[i - 1];
233  y1 = y2;
234  y2 = y3;
235  }
236  }
237  thinned->points[j++] = ptwXY1->points[length - 1];
238 
239  if( ptwXY1->interpolation != ptwXY_interpolationFlat ) { /* Now call ptwXY_thin2 for more thinning. */
240  length = thinned->length = j;
241  if( ( thin = (char *) nfu_calloc( 1, (size_t) length ) ) == NULL ) goto Err;
242  if( ( *status = ptwXY_thin2( thinned, thin, accuracy, 0, length - 1 ) ) != nfu_Okay ) goto Err;
243  for( j = 1; j < length; j++ ) if( thin[j] != 0 ) break;
244  for( i = j + 1; i < length; i++ ) {
245  if( thin[i] == 0 ) {
246  thinned->points[j] = thinned->points[i];
247  j++;
248  }
249  }
250  nfu_free( thin );
251  }
252  thinned->length = j;
253 
254  return( thinned );
255 
256 Err:
257  ptwXY_free( thinned );
258  if( thin != NULL ) nfu_free( thin );
259  return( NULL );
260 }
261 /*
262 ************************************************************
263 */
264 static nfu_status ptwXY_thin2( ptwXYPoints *thinned, char *thin, double accuracy, int64_t i1, int64_t i2 ) {
265 
266  int64_t i, iMax = 0;
267  double y, s, dY, dYMax = 0., dYR, dYRMax = 0;
268  double x1 = thinned->points[i1].x, y1 = thinned->points[i1].y, x2 = thinned->points[i2].x, y2 = thinned->points[i2].y;
269  nfu_status status;
270 
271  if( i1 + 1 >= i2 ) return( nfu_Okay );
272  for( i = i1 + 1; i < i2; i++ ) {
273  if( ( status = ptwXY_interpolatePoint( thinned->interpolation, thinned->points[i].x, &y, x1, y1, x2, y2 ) ) != nfu_Okay ) return( status );
274  s = 0.5 * ( std::fabs( y ) + std::fabs( thinned->points[i].y ) );
275  dY = std::fabs( y - thinned->points[i].y );
276  dYR = 0;
277  if( s != 0 ) dYR = dY / s;
278  if( ( dYR > dYRMax ) || ( ( dYR >= 0.9999 * dYRMax ) && ( dY > dYMax ) ) ) { /* The choice of 0.9999 is not exact science. */
279  iMax = i;
280  if( dY > dYMax ) dYMax = dY;
281  if( dYR > dYRMax ) dYRMax = dYR;
282  }
283  }
284  if( dYRMax < accuracy ) {
285  for( i = i1 + 1; i < i2; i++ ) thin[i] = 1; }
286  else {
287  if( ( status = ptwXY_thin2( thinned, thin, accuracy, i1, iMax ) ) != nfu_Okay ) return( status );
288  status = ptwXY_thin2( thinned, thin, accuracy, iMax, i2 );
289  }
290  return( status );
291 }
292 /*
293 ************************************************************
294 */
295 static double ptwXY_thicken_linear_dx( int sectionSubdivideMax, double dxMax, double x1, double x2 ) {
296 
297  int ndx;
298  double dx = x2 - x1, dndx;
299 
300  if( dxMax == 0. ) {
301  dx = ( x2 - x1 ) / sectionSubdivideMax; }
302  else {
303  dndx = dx / dxMax;
304  ndx = (int) dndx;
305  if( ( dndx - ndx ) > 1e-6 ) ndx++;
306  if( ndx > sectionSubdivideMax ) ndx = sectionSubdivideMax;
307  if( ndx > 0 ) dx /= ndx;
308  }
309 
310  return( dx );
311 }
312 /*
313 ************************************************************
314 */
316 /*
317 c Remove extra zeros at beginning and end.
318 */
319 
320  int64_t i, i1, i2;
321  nfu_status status;
322 
323  if( ptwXY->status != nfu_Okay ) return( ptwXY->status );
324  if( ( status = ptwXY_simpleCoalescePoints( ptwXY ) ) != nfu_Okay ) return( status );
325  for( i1 = 0; i1 < ptwXY->length; i1++ ) {
326  if( ptwXY->points[i1].y != 0 ) break;
327  }
328  if( i1 > 0 ) i1--;
329  for( i2 = ptwXY->length - 1; i2 >= 0; i2-- ) {
330  if( ptwXY->points[i2].y != 0 ) break;
331  }
332  i2++;
333  if( i2 < ptwXY->length ) i2++;
334  if( i2 > i1 ) {
335  if( i1 > 0 ) {
336  for( i = i1; i < i2; i++ ) ptwXY->points[i - i1] = ptwXY->points[i];
337  }
338  ptwXY->length = i2 - i1; }
339  else if( i2 < i1 ) { /* Remove all zeros between endpoints. */
340  ptwXY->points[1] = ptwXY->points[ptwXY->length - 1];
341  ptwXY->length = 2;
342  }
343 
344  return( nfu_Okay );
345 }
346 /*
347 ************************************************************
348 */
349 ptwXYPoints *ptwXY_union( ptwXYPoints *ptwXY1, ptwXYPoints *ptwXY2, nfu_status *status, int unionOptions ) {
350 
351  int64_t overflowSize, i, i1 = 0, i2 = 0, n1 = ptwXY1->length, n2 = ptwXY2->length, length;
352  int fillWithFirst = unionOptions & ptwXY_union_fill, trim = unionOptions & ptwXY_union_trim;
353  ptwXYPoints *n;
354  double x1 = 0., x2 = 0., y1 = 0., y2 = 0., y, biSectionMax, accuracy;
355 
356  if( ( *status = ptwXY1->status ) != nfu_Okay ) return( NULL );
357  if( ( *status = ptwXY2->status ) != nfu_Okay ) return( NULL );
358  *status = nfu_otherInterpolation;
359  if( ptwXY1->interpolation == ptwXY_interpolationOther ) return( NULL );
360 /*
361 * Many other routines use the fact that ptwXY_union calls ptwXY_coalescePoints for ptwXY1 and ptwXY2 so do not change it.
362 */
363  if( ( *status = ptwXY_simpleCoalescePoints( ptwXY1 ) ) != nfu_Okay ) return( NULL );
364  if( ( *status = ptwXY_simpleCoalescePoints( ptwXY2 ) ) != nfu_Okay ) return( NULL );
365 
366  if( ( n1 == 1 ) || ( n2 == 1 ) ) {
367  *status = nfu_tooFewPoints;
368  return( NULL );
369  }
370  if( trim ) {
371  if( n1 > 0 ) {
372  if( n2 > 0 ) {
373  if( ptwXY1->points[0].x < ptwXY2->points[0].x ) {
374  while( i1 < n1 ) { // Loop checking, 11.05.2015, T. Koi
375  if( ptwXY1->points[i1].x >= ptwXY2->points[0].x ) break;
376  if( fillWithFirst ) {
377  if( i1 < ( ptwXY1->length - 1 ) ) {
378  x1 = ptwXY1->points[i1].x;
379  y1 = ptwXY1->points[i1].y;
380  x2 = ptwXY1->points[i1+1].x;
381  y2 = ptwXY1->points[i1+1].y;
382  }
383  }
384  i1++;
385  } }
386  else {
387  while( i2 < n2 ) { // Loop checking, 11.06.2015, T. Koi
388  if( ptwXY2->points[i2].x >= ptwXY1->points[0].x ) break;
389  i2++;
390  }
391  }
392  if( ptwXY1->points[n1-1].x > ptwXY2->points[n2-1].x ) {
393  while( i1 < n1 ) { // Loop checking, 11.06.2015, T. Koi
394  if( ptwXY1->points[n1-1].x <= ptwXY2->points[n2-1].x ) break;
395  n1--;
396  } }
397  else {
398  while( i2 < n2 ) { // Loop checking, 11.06.2015, T. Koi
399  if( ptwXY2->points[n2-1].x <= ptwXY1->points[n1-1].x ) break;
400  n2--;
401  }
402  } }
403  else {
404  n1 = 0;
405  } }
406  else {
407  n2 = 0;
408  }
409  }
410  overflowSize = ptwXY1->overflowAllocatedSize;
411  if( overflowSize < ptwXY2->overflowAllocatedSize ) overflowSize = ptwXY2->overflowAllocatedSize;
412  length = ( n1 - i1 ) + ( n2 - i2 );
413  if( length == 0 ) length = ptwXY_minimumSize;
414  biSectionMax = ptwXY1->biSectionMax;
415  if( biSectionMax < ptwXY2->biSectionMax ) biSectionMax = ptwXY2->biSectionMax;
416  accuracy = ptwXY1->accuracy;
417  if( accuracy < ptwXY2->accuracy ) accuracy = ptwXY2->accuracy;
418  n = ptwXY_new( ptwXY1->interpolation, NULL, biSectionMax, accuracy, length, overflowSize, status, ptwXY1->userFlag );
419  if( n == NULL ) return( NULL );
420 
421  for( i = 0; ( i1 < n1 ) && ( i2 < n2 ); i++ ) {
422  y = 0.;
423  if( ptwXY1->points[i1].x <= ptwXY2->points[i2].x ) {
424  n->points[i].x = ptwXY1->points[i1].x;
425  if( fillWithFirst ) {
426  y = ptwXY1->points[i1].y;
427  if( i1 < ( ptwXY1->length - 1 ) ) {
428  x1 = ptwXY1->points[i1].x;
429  y1 = ptwXY1->points[i1].y;
430  x2 = ptwXY1->points[i1+1].x;
431  y2 = ptwXY1->points[i1+1].y; }
432  else {
433  y1 = 0.;
434  y2 = 0.;
435  }
436  }
437  if( ptwXY1->points[i1].x == ptwXY2->points[i2].x ) i2++;
438  i1++; }
439  else {
440  n->points[i].x = ptwXY2->points[i2].x;
441  if( fillWithFirst && ( ( y1 != 0. ) || ( y2 != 0. ) ) ) {
442  if( ( *status = ptwXY_interpolatePoint( ptwXY1->interpolation, ptwXY2->points[i2].x, &y, x1, y1, x2, y2 ) ) != nfu_Okay ) {
443  ptwXY_free( n );
444  return( NULL );
445  }
446  }
447  i2++;
448  }
449  n->points[i].y = y;
450  }
451 
452  y = 0.;
453  for( ; i1 < n1; i1++, i++ ) {
454  n->points[i].x = ptwXY1->points[i1].x;
455  if( fillWithFirst ) y = ptwXY1->points[i1].y;
456  n->points[i].y = y;
457  }
458  for( ; i2 < n2; i2++, i++ ) {
459  n->points[i].x = ptwXY2->points[i2].x;
460  if( fillWithFirst && trim && ( n->points[i].x <= x2 ) ) {
461  if( ( *status = ptwXY_interpolatePoint( ptwXY1->interpolation, n->points[i].x, &y, x1, y1, x2, y2 ) ) != nfu_Okay ) {
462  ptwXY_free( n );
463  return( NULL );
464  }
465  }
466  n->points[i].y = y;
467  }
468  n->length = i;
469 
470  if( unionOptions & ptwXY_union_mergeClosePoints ) {
471  if( ( *status = ptwXY_mergeClosePoints( n, 4 * DBL_EPSILON ) ) != nfu_Okay ) {
472  ptwXY_free( n );
473  return( NULL );
474  }
475  }
476  return( n );
477 }
478 /*
479 ************************************************************
480 */
481 nfu_status ptwXY_scaleOffsetXAndY( ptwXYPoints *ptwXY, double xScale, double xOffset, double yScale, double yOffset ) {
482 
483  int64_t i1, length = ptwXY->length;
484  ptwXYPoint *p1;
485  nfu_status status;
486 
487  if( ptwXY->status != nfu_Okay ) return( ptwXY->status );
488  if( xScale == 0 ) return( nfu_XNotAscending );
489 
490  if( ( status = ptwXY_simpleCoalescePoints( ptwXY ) ) != nfu_Okay ) return( status );
491 
492  for( i1 = 0, p1 = ptwXY->points; i1 < length; i1++, p1++ ) {
493  p1->x = xScale * p1->x + xOffset;
494  p1->y = yScale * p1->y + yOffset;
495  }
496 
497  if( xScale < 0 ) {
498  int64_t length_2 = length / 2;
499  ptwXYPoint tmp, *p2;
500 
501  for( i1 = 0, p1 = ptwXY->points, p2 = &(ptwXY->points[length-1]); i1 < length_2; i1++ ) {
502  tmp = *p1;
503  *p1 = *p2;
504  *p2 = tmp;
505  }
506  }
507 
508  return( ptwXY->status );
509 }
510 
511 #if defined __cplusplus
512 }
513 #endif