~ [ source navigation ] ~ [ diff markup ] ~ [ identifier search ] ~ [ freetext search ] ~ [ file search ] ~

Linux Cross Reference
Tina5/tina-libs/tina/math/mathUtil_hist.c

Version: ~
Architecture: ~ [ i386 ] ~ [ alpha ] ~ [ m68k ] ~ [ mips ] ~ [ ppc ] ~ [ sparc ] ~ [ sparc64 ] ~

  1 /**********
  2  *
  3  * Copyright (c) 2003, Division of Imaging Science and Biomedical Engineering,
  4  * University of Manchester, UK.  All rights reserved.
  5  * 
  6  * Redistribution and use in source and binary forms, with or without modification, 
  7  * are permitted provided that the following conditions are met:
  8  * 
  9  *   . Redistributions of source code must retain the above copyright notice, 
 10  *     this list of conditions and the following disclaimer.
 11  *    
 12  *   . Redistributions in binary form must reproduce the above copyright notice,
 13  *     this list of conditions and the following disclaimer in the documentation 
 14  *     and/or other materials provided with the distribution.
 15  * 
 16  *   . Neither the name of the University of Manchester nor the names of its
 17  *     contributors may be used to endorse or promote products derived from this 
 18  *     software without specific prior written permission.
 19  * 
 20  * 
 21  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" 
 22  * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 
 23  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 
 24  * ARE DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE 
 25  * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR 
 26  * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF 
 27  * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
 28  * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN 
 29  * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) 
 30  * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE 
 31  * POSSIBILITY OF SUCH DAMAGE.
 32  *
 33  **********
 34  *
 35  * Program :    TINA
 36  * File    :  $Source: /home/tina/cvs/tina-libs/tina/math/mathUtil_hist.c,v $
 37  * Date    :  $Date: 2008/10/02 11:53:05 $
 38  * Version :  $Revision: 1.15 $
 39  * CVS Id  :  $Id: mathUtil_hist.c,v 1.15 2008/10/02 11:53:05 neil Exp $
 40  *
 41  * Author  :  Legacy TINA
 42  *
 43  * Notes :
 44  *
 45  *********
 46 */
 47 /** 
 48  *  @file       
 49  *  @brief   Histogram creation and utility functions.
 50  *
 51  *  Functions include;
 52  * 
 53  *  - creation and freeing of 1D and 2D histograms.
 54  *  - Histogram filling.
 55  *  - Histogram formatting/printing.
 56  *  - Calculation of various statistical properties of the histogram. 
 57 */
 58 
 59 #include "mathUtil_hist.h"
 60 
 61 #if HAVE_CONFIG_H
 62 #include <config.h>
 63 #endif
 64 
 65 #include <stdio.h>
 66 #include <stdlib.h>
 67 #include <math.h>
 68 #include <float.h>
 69 #include <tina/sys/sysDef.h>
 70 #include <tina/sys/sysPro.h>
 71 #include <tina/math/mathDef.h>
 72 #include <tina/math/mathPro.h>
 73 #include <tina/math/math_UtilDef.h>
 74 
 75 double          (*mfunc) (int order, double *params, float x) = NULL;
 76 double          (*merror_func) (shistogram * ph, float x) = NULL;
 77 
 78 static void     herror(char *str);
 79 
 80 static shistogram **hist=NULL;                  /* static data! */
 81 static shistogram *mph;                         /* static data! */
 82 
 83 
 84 shistogram    **hist_vec()
 85 {
 86         int i;
 87         if (hist == NULL)
 88         {
 89             hist = (shistogram **) ralloc(NHIST * sizeof(void *));
 90             for (i=0; i< NHIST ;i++)
 91             {
 92                  hist[i] = NULL;
 93             }
 94         }
 95 
 96         return (hist);
 97 }
 98 
 99 
100 static double **maketarray(int l_array, int w_array, float gl)
101 {
102         double        **array;
103         short           i, j;
104 
105         if ((array = (double **) ralloc((unsigned) l_array * sizeof(double *)))
106             == NULL)
107                 return (NULL);
108         if ((array[0] = (double *) ralloc((unsigned) l_array * w_array * sizeof(double)))
109             == NULL)
110                 return (void *) (NULL);
111         for (i = 1; i < l_array; ++i)
112                 array[i] = array[i - 1] + w_array;
113         for (i = 0; i < l_array; ++i)
114                 for (j = 0; j < w_array; j++)
115                         array[i][j] = gl;
116         return (array);
117 }
118 
119 static void     freetarray(double **array)
120 {
121         if (array == NULL)
122                 return;
123         if (array[0] != NULL)
124                 rfree((void *) array[0]);
125         if (array != NULL)
126                 rfree((void *) array);
127 }
128 
129 static float  **makeharray(int l_array, int w_array, float gl)
130 {
131         float         **array;
132         short           i, j;
133         if ((array = (float **) ralloc(l_array * sizeof(float *))) == NULL)
134                 herror("no room for array");
135         if ((array[0] = (float *) ralloc(l_array * w_array * sizeof(float))) == NULL)
136                 herror("no room for array");
137         for (i = 1; i < l_array; ++i)
138                 array[i] = array[i - 1] + w_array;
139         for (i = 0; i < l_array; ++i)
140                 for (j = 0; j < w_array; j++)
141                         array[i][j] = gl;
142         return (array);
143 }
144 
145 static void     freeharray(float **array)
146 {
147         if (array == NULL)
148                 return;
149         if (array[0] != NULL)
150                 rfree(array[0]);
151         if (array != NULL)
152                 rfree(array);
153 }
154 
155 shistogram     *hbook1(int id, char *title, float xmin, float xmax, int xbins)
156 {
157         shistogram     *ph;
158 
159         if (id >= NHIST)
160         {
161                 herror("cannot book histogram \n");
162                 return (NULL);
163         }
164         ph = (shistogram *) ralloc(sizeof(shistogram));
165         ph->id = id;
166         ph->title = title;
167         ph->type = 1.0;
168         ph->shFuncSuper = NULL;
169         ph->xmax = xmax;
170         ph->xmin = xmin;
171         ph->contents = 0.0;
172         ph->mean = 0.0;
173         ph->mean2 = 0.0;
174         ph->under = 0.0;
175         ph->over = 0.0;
176         ph->entries = 0.0;
177         ph->xbins = xbins;
178         ph->ybins = 1;
179         ph->xincr = (xmax - xmin) / (float) xbins;
180         ph->array = makeharray(1, xbins, 0.0);
181         return (ph);
182 }
183 
184 shistogram     *hfree(shistogram * ph)
185 {
186         if (ph != NULL)
187         {
188                 freeharray(ph->array);
189                 if (ph->par != NULL)
190                         rfree(ph->par);
191                 rfree(ph);
192         }
193         return (NULL);
194 }
195 
196 shistogram     *hbook2(int id, char *title, float xmin, float xmax, int xbins,
197                                        float ymin, float ymax, int ybins)
198 {
199         shistogram     *ph;
200 
201         if (id >= NHIST)
202         {
203                 herror("cannot book histogram \n");
204                 return (NULL);
205         }
206         ph = (shistogram *) ralloc(sizeof(shistogram));
207         ph->id = id;
208         ph->title = title;
209         ph->type = 2;
210         ph->shFuncSuper = NULL;
211         ph->xmax = xmax;
212         ph->xmin = xmin;
213         ph->ymin = ymin;
214         ph->ymax = ymax;
215         ph->contents = 0.0;
216         ph->entries = 0.0;
217         ph->mean = 0.0;
218         ph->under = 0.0;
219         ph->over = 0.0;
220         ph->above = 0.0;
221         ph->below = 0.0;
222         ph->xbins = xbins;
223         ph->ybins = ybins;
224         ph->xincr = (xmax - xmin) / (float) xbins;
225         ph->yincr = (ymax - ymin) / (float) ybins;
226         ph->array = makeharray(ybins, xbins, 0.0);
227         return (ph);
228 }
229 
230 float           hfill1(shistogram * ph, float x, float w)
231 {
232         int             i;
233         if (ph == NULL)
234                 return (0.0);
235         i = floor((x - ph->xmin) / ph->xincr);
236         if (x < ph->xmax && x >= ph->xmin)
237         {
238                 if (w == 0.0)
239                         return (ph->array[0][i]);
240                 ph->contents += w;
241                 ph->entries++;
242                 ph->mean += w * (x - ph->mean) / ph->contents;
243                 ph->mean2 += w * (x * x - ph->mean2) / ph->contents;
244                 /* correct numerical instability on first time through */
245                 if (ph->mean * ph->mean > ph->mean2)
246                         ph->mean2 = ph->mean * ph->mean;
247                 return (ph->array[0][i] += w);
248         } else
249         {
250                 if (x > ph->xmax)
251                         ph->over += w;
252                 else
253                         ph->under += w;
254                 return (0.0);
255         }
256 }
257 
258 float           hfill1s(shistogram * ph, float x, float w)
259 {
260         int             i;
261         double          fac;
262         if (ph == NULL)
263                 return (0.0);
264         i = floor((x - ph->xmin) / ph->xincr + 0.5);
265         if (x < ph->xmax && x >= ph->xmin)
266         {
267                 if (w == 0.0)
268                         return (ph->array[0][i]);
269                 ph->contents += w;
270                 ph->entries++;
271                 ph->mean += w * (x - ph->mean) / ph->contents;
272                 ph->mean2 += w * (x * x - ph->mean2) / ph->contents;
273                 /* correct numerical instability on first time through */
274                 if (ph->mean * ph->mean > ph->mean2)
275                         ph->mean2 = ph->mean * ph->mean;
276                 fac = (x - i * ph->xincr - ph->xmin) / ph->xincr;
277                 if (i > 0)
278                         ph->array[0][i - 1] += w * (0.5 - fac);
279                 if (i < ph->xbins)
280                         ph->array[0][i] += w * (0.5 + fac);
281                 return (ph->array[0][i]);
282         } else
283         {
284                 if (x > ph->xmax)
285                         ph->over += w;
286                 else
287                         ph->under += w;
288                 return (0.0);
289         }
290 }
291 
292 float           hfill2(shistogram * ph, float x, float y, float w)
293 {
294         int             i, j;
295         if (ph == NULL || w == 0.0)
296                 return (0.0);
297         i = floor((x - ph->xmin) / ph->xincr);
298         j = floor((y - ph->ymin) / ph->yincr);
299         if (x < ph->xmax && x >= ph->xmin && y < ph->ymax && y >= ph->ymin)
300         {
301                 ph->contents += w;
302                 ph->entries++;
303                 ph->mean += w * (x - ph->mean) / ph->contents;
304                 return (ph->array[j][i] += w);
305         } else
306         {
307                 if (x < ph->xmin)
308                         return (ph->under += w);
309                 if (x >= ph->xmax)
310                         return (ph->over += w);
311                 if (y < ph->ymin)
312                         return (ph->below += w);
313                 if (y >= ph->ymax)
314                         return (ph->above += w);
315         }
316         return (0.0);           /* should never be executed */
317 }
318 
319 
320 float hfill2s(shistogram *ph,float x,float y,float w)
321 {
322    /* Histogram filling function written specifically for the mutual information 
323       coregistration: subpixel interpolation (linear) on input and (quadratic) 
324       on output */
325 
326    int i,j, n, m;
327    float xs, ys, a, b, c, d, e, f;
328    double ifac, jfac, inter, pixval[3][3];
329    if (ph==NULL) return(0.0);
330    i = floor((x-ph->xmin)/ph->xincr + 0.5);
331    j = floor((y-ph->ymin)/ph->yincr + 0.5);
332    if (x<ph->xmax && x>=ph->xmin && y<ph->ymax && y>=ph->ymin)
333    {
334         ph->contents+=w;
335         if(w!=0.0) ph->entries++;
336         ph->mean += w*(x-ph->mean)/ph->contents;
337         ifac = (x-i*ph->xincr-ph->xmin)/ph->xincr;
338         jfac = (y-j*ph->yincr-ph->ymin)/ph->yincr;
339 
340         /* PAB updates 30/5/2003: make sure that all entries go into 
341          the histogram somewhere */
342 
343         if (i>0&&j>0) 
344         {
345                 ph->array[j-1][i-1]+=w*(0.5-ifac)*(0.5-jfac); 
346         }
347         else
348         {
349                 if(j==0)
350                 {
351                         ph->array[j][i-1]+=w*(0.5-ifac)*(0.5-jfac); 
352                 }
353                 else
354                 {
355                         ph->array[j-1][i]+=w*(0.5-ifac)*(0.5-jfac);     
356                 }
357         }
358 
359         if (i<ph->xbins&&j>0) 
360         {
361                 ph->array[j-1][i]+=w*(0.5-jfac)*(0.5+ifac); 
362         }
363         else
364         {
365                 if(j==0)
366                 {
367                         ph->array[j][i]+=w*(0.5-jfac)*(0.5+ifac); 
368                 }
369                 else
370                 {
371                         ph->array[j-1][i-1]+=w*(0.5-jfac)*(0.5+ifac); 
372                 }
373         }
374 
375         if (i>0&&j<ph->ybins)
376         {
377                 ph->array[j][i-1]+=w*(0.5+jfac)*(0.5-ifac);
378         }
379         else
380         {
381                 if(i==0)
382                 {
383                         ph->array[j][i]+=w*(0.5+jfac)*(0.5-ifac);
384                 }
385                 else
386                 {
387                         ph->array[j-1][i-1]+=w*(0.5+jfac)*(0.5-ifac);
388                 }
389         }
390 
391         if (i<ph->xbins&&j<ph->ybins)
392         {
393                 ph->array[j][i]+=w*(0.5+ifac)*(0.5+jfac);
394         }
395         else
396         {
397                 if(i==ph->xbins)
398                 {
399                         ph->array[j][i-1]+=w*(0.5+ifac)*(0.5+jfac);
400                 }
401                 else
402                 {
403                         ph->array[j-1][i]+=w*(0.5+ifac)*(0.5+jfac);
404                 }
405         }
406 
407    /* Subpixel interpolate a return value, using a copy of the im_sub_pixqf code
408       from imget.c */
409 
410    i = floor((x-ph->xmin)/ph->xincr)-1;
411    j = floor((y-ph->ymin)/ph->yincr)-1;
412 
413     xs = (float)((x-(i+1.5)*ph->xincr-ph->xmin)/ph->xincr);
414     ys = (float)((y-(j+1.5)*ph->yincr-ph->ymin)/ph->yincr);
415 
416     if (j < 0|| j > ph->ybins - 3
417         || i < 0 || i >  ph->xbins - 3)
418         return (0);
419     for (n = 0; n < 3; n++)
420     {
421         for (m = 0; m < 3; m++)
422         {
423             pixval[n][m] = (double)(ph->array[(j+n)][(i+m)]);
424         }
425     }
426 
427     a = (float)pixval[1][1];
428     b = (float)((pixval[0][2] - pixval[0][0]
429         + pixval[1][2] - pixval[1][0]
430         + pixval[2][2] - pixval[2][0]) / 6.0);
431     c = (float)((pixval[2][0] - pixval[0][0]
432         + pixval[2][1] - pixval[0][1]
433         + pixval[2][2] - pixval[0][2]) / 6.0);
434     d = (float)((pixval[0][0] - 2.0 * pixval[0][1] + pixval[0][2]
435         + 3.0 * pixval[1][0] - 6.0 * pixval[1][1] + 3.0 * pixval[1][2]
436         + pixval[2][0] - 2.0 * pixval[2][1] + pixval[2][2]) / 10.0);
437     e = (float)((pixval[0][0] - pixval[2][0]
438         + pixval[2][2] - pixval[0][2]) / 4.0);
439     f = (float)((pixval[0][0] + 3.0 * pixval[0][1] + pixval[0][2]
440         - 2.0 * pixval[1][0] - 6.0 * pixval[1][1] - 2.0 * pixval[1][2]
441         + pixval[2][0] + 3.0 * pixval[2][1] + pixval[2][2]) / 10.0);
442 
443     xs *=(float)2.0;
444     ys *=(float)2.0;
445     inter = a + 0.5*((f*ys+e*xs+c)*ys + (d*xs+b)*xs);
446 
447     return (inter);
448    }
449    else 
450    {
451       if(x<ph->xmin) return(ph->under +=w);
452       if(x>=ph->xmax) return(ph->over +=w);
453       if(y<ph->ymin) return(ph->below += w);
454       if(y>=ph->ymax) return(ph->above += w);
455       return 0.0;
456    }
457 }
458 
459 
460 void            hprint(FILE * fp, shistogram * ph)
461 {
462         int             i, j;
463         float           scale = 1.0, max, variance, y;
464         char            temp[20];
465         if (ph == NULL)
466         {
467                 fprintf(fp, "histogram empty! \n");
468                 return;
469         }
470         if (ph->entries == 0)
471                 fprintf(fp, " h id %d %s empty \n\n", ph->id, ph->title);
472         else
473                 fprintf(fp, " h id %d %s\n\n", ph->id, ph->title);
474 
475         if (ph->type == 1)
476         {
477                 if (ph->entries != 0)
478                 {
479                         max = 0.0;
480                         for (i = 0; i < ph->xbins; i++)
481                         {
482                                 if (ph->array[0][i] > max)
483                                         max = ph->array[0][i];
484                         }
485                         scale = max / 100.0;
486                         scale = pow(10.0, (double) ((int) (log10((double) scale))));
487                         if (max > 50.0 * scale)
488                                 scale *= 2.0;
489                         if (max > 50.0 * scale)
490                                 scale *= 2.5;
491                         if (max > 50.0 * scale)
492                                 scale *= 2.0;
493                         if (max > 50.0 * scale)
494                                 scale *= 2.0;
495 
496                         for (i = 0; i < ph->xbins; i++)
497                         {
498                                 fprintf(fp, " %6.2f ", ph->xmin + i * ph->xincr);
499                                 if (ph->shFuncSuper == NULL)
500                                 {
501                                         for (j = 0; (float) j < ph->array[0][i] / scale; j++)
502                                         {
503                                                 fprintf(fp, "x");
504                                         }
505                                 } else
506                                 {
507                                         y = ph->shFuncSuper(ph->npar, ph->par, (ph->xmin + (i + 0.5) * ph->xincr));
508                                         if (y / scale > 60.0)
509                                                 y = scale * 60.0;
510 
511                                         for (j = 0; (float) j < ph->array[0][i] / scale
512                                              || (float) j <= y / scale; j++)
513                                         {
514                                                 if ((int) (y / scale) == j)
515                                                         fprintf(fp, "*");
516                                                 else if ((float) j + 0.5 < ph->array[0][i] / scale)
517                                                         fprintf(fp, "x");
518                                                 else
519                                                         fprintf(fp, " ");
520                                         }
521                                 }
522                                 fprintf(fp, "\n");
523                         }
524                         fprintf(fp, "\n");
525                 }
526                 fprintf(fp, " h id %d entries %d contents %5.3f underflows %5.3f overflows %5.3f\n",
527                     ph->id, ph->entries, ph->contents, ph->under, ph->over);
528                 variance = (float) sqrt(fabs(ph->mean2 - ph->mean * ph->mean));
529                 fprintf(fp, " h id %d mean  %5.3f var(**0.5) %5.3f scale %5.3f res %5.3f run %5.3f grad %5.3f \n",
530                         ph->id, ph->mean, variance, scale, hlnorm(ph, 2), hrunstat(ph), hgradstat(ph));
531                 fprintf(fp, "\n");
532                 fprintf(fp, "\n");
533         } else if (ph->type == 2)
534         {
535                 if (ph->entries != 0)
536                 {
537                         for (j = ph->ybins - 1; j >= 0; j--)
538                         {
539                                 fprintf(fp, " %6.2f ", ph->ymin + j * ph->yincr);
540                                 for (i = 0; i < ph->xbins; i++)
541                                 {
542                                         if (ph->array[j][i] >= 1.0 && ph->array[j][i] < 10.0)
543                                                 fprintf(fp, "%d", (int) ph->array[j][i]);
544                                         else if (ph->array[j][i] == 0.0)
545                                                 fprintf(fp, " ");
546                                         else if (ph->array[j][i] < 1.0 && ph->array[j][i] > 0.0)
547                                                 fprintf(fp, ".");
548                                         else if (ph->array[j][i] >= 200.0)
549                                                 fprintf(fp, "*");
550                                         else if (ph->array[j][i] >= 150.0)
551                                                 fprintf(fp, "H");
552                                         else if (ph->array[j][i] >= 100.0)
553                                                 fprintf(fp, "G");
554                                         else if (ph->array[j][i] >= 65.0)
555                                                 fprintf(fp, "F");
556                                         else if (ph->array[j][i] >= 45.0)
557                                                 fprintf(fp, "E");
558                                         else if (ph->array[j][i] >= 30.0)
559                                                 fprintf(fp, "D");
560                                         else if (ph->array[j][i] >= 20.0)
561                                                 fprintf(fp, "C");
562                                         else if (ph->array[j][i] >= 15.0)
563                                                 fprintf(fp, "B");
564                                         else if (ph->array[j][i] >= 10.0)
565                                                 fprintf(fp, "A");
566                                         else if (ph->array[j][i] < 0.0)
567                                                 fprintf(fp, "-");
568                                         else
569                                                 fprintf(fp, "?");
570                                 }
571                                 fprintf(fp, "\n");
572                         }
573                         sprintf(temp, "%4.2f", ph->xmin);
574                         for (j = 0; temp[j] != '\0'; j++)
575                         {
576                                 fprintf(fp, "        ");
577                                 for (i = 0; i < ph->xbins; i++)
578                                 {
579                                         sprintf(temp, "%4.2f", ph->xmin + i * ph->xincr);
580                                         fprintf(fp, "%c", temp[j]);
581                                 }
582                                 fprintf(fp, "\n");
583                         }
584                         fprintf(fp, "\n");
585                         fprintf(fp, " Key: .<1.0,  A<15.0,  B<20.0,  C<30.0,  D<45.0,\n");
586                         fprintf(fp, "      E<65.0  F<100.0, G<150.0, H<200.0, *<Infty.\n");
587                 }
588                 fprintf(fp, " entries %d contents %6.2f\n", ph->entries, ph->contents);
589                 fprintf(fp, " under%6.2f over %6.2f above %6.2f below %6.2f\n",
590                         ph->under, ph->over, ph->above, ph->below);
591                 fprintf(fp, "\n");
592         }
593 }
594 
595 double          hresidual(shistogram * ph, int i)
596 {
597         double          y, residual;
598         y = ph->shFuncSuper(ph->npar, ph->par, (ph->xmin + (i + 0.5) * ph->xincr));
599         residual = y - ph->array[0][i];
600         return (residual);
601 }
602 
603 double          hlnorm(shistogram * ph, int norm)
604 {
605         int             i;
606         double          residual = 0.0;
607         if (ph->shFuncSuper == NULL)
608                 return (0.0);
609         for (i = 1; i < ph->xbins - 1; i++)
610         {
611                 if (norm == 1)
612                         residual += fabs(hresidual(ph, i));
613                 else if (norm == 2)
614                         residual += hresidual(ph, i) * hresidual(ph, i);
615         }
616         residual /= (float) (ph->xbins - 2);
617 
618         if (norm == 1)
619                 return (residual);
620         else
621                 return (sqrt(residual));
622 }
623 
624 double          hgradstat(shistogram * ph)
625 {
626         int             i;
627         double          y, old_y, data, old_data;
628         double          gradtest;
629         int             count = 0;
630         if (ph->shFuncSuper == NULL)
631                 return (0.0);
632         old_y = ph->shFuncSuper(ph->npar, ph->par, (ph->xmin + (0.5) * ph->xincr));
633         old_data = ph->array[0][0];
634 
635         for (i = 1; i < ph->xbins; i++)
636         {
637                 y = ph->shFuncSuper(ph->npar, ph->par, (ph->xmin + (i + 0.5) * ph->xincr));
638                 data = ph->array[0][i];
639                 gradtest = (y - old_y) * (data - old_data);
640                 if (gradtest > 0.0)
641                 {
642                         count++;
643                 }
644                 old_y = y;
645                 old_data = data;
646         }
647 
648         return ((double) count / ((float) ph->xbins - 1));
649 }
650 
651 double          hrunstat(shistogram * ph)
652 {
653         int             i, count[8] = {0, 0, 0, 0, 0, 0, 0, 0};
654         int             ncount = 0;
655         double          new_residual, residual = 0.0, y;
656         double          bhatt = 0.0, temp = 1.0;
657 
658         if (ph->shFuncSuper == NULL)
659                 return (0.0);
660         /* analyse run count */
661         for (i = 0; i < ph->xbins; i++)
662         {
663                 y = ph->shFuncSuper(ph->npar, ph->par, (ph->xmin + (i + 0.5) * ph->xincr));
664                 new_residual = y - ph->array[0][i];
665                 if ((residual < 0.0 && new_residual < 0.0)
666                     || (residual >= 0.0 && new_residual >= 0.0))
667                 {
668                         if (i != 0)
669                                 ncount++;
670                 } else
671                 {
672                         ncount = 0;
673                 }
674                 if (ncount > 7)
675                         ncount = 7;
676                 count[ncount]++;
677                 residual = new_residual;
678         }
679         /* compare frequency distribution against expected using B statistic */
680         for (i = 0; i < 7; i++)
681         {
682                 temp *= 0.5;
683                 bhatt += sqrt(temp) * sqrt((float) count[i]);
684         }
685         /* finish off last bin */
686         bhatt += sqrt(temp) * sqrt((float) count[7]);
687         return (bhatt / sqrt((float) ph->xbins));
688 }
689 
690 void            histdo(FILE * fp)
691 {
692         int             i;
693         if (fp == NULL)
694                 return;
695         for (i = 0; i < NHIST; i++)
696         {
697                 if (hist[i] != NULL)
698                         hprint(fp, hist[i]);
699         }
700 }
701 
702 void            hreset(shistogram * ph)
703 {
704         int             i, j;
705 
706         if (ph == NULL)
707         {
708                 for (i = 0; i < NHIST; i++)
709                 {
710                         if (hist[i] != NULL)
711                                 hreset(hist[i]);
712                 }
713         } else
714         {
715                 for (j = 0; j < ph->ybins; j++)
716                 {
717                         for (i = 0; i < ph->xbins; i++)
718                         {
719                                 ph->array[j][i] = 0.0;
720                         }
721                 }
722                 ph->contents = 0.0;
723                 ph->entries = 0.0;
724                 ph->under = 0.0;
725                 ph->over = 0.0;
726                 ph->above = 0.0;
727                 ph->below = 0.0;
728                 ph->mean = 0.0;
729                 ph->mean2 = 0.0;
730         }
731 }
732 
733 void            hpxprint(FILE * fp, shistogram * ph)
734 {
735         int             i, j;
736         char           *title;
737         shistogram     *temp = NULL;
738 
739         if (ph == NULL)
740         {
741                 return;
742         } else
743         {
744                 title = (char *) ralloc(256 * sizeof(char));
745                 sprintf(title, " X projection > %s ", ph->title);
746                 temp = hbook1(NHIST - 1, title, ph->xmin, ph->xmax, ph->xbins);
747                 for (j = 0; j < ph->ybins; j++)
748                 {
749                         for (i = 0; i < ph->xbins; i++)
750                         {
751                                 hfill1(temp, ph->xmin + (i + 0.5) * ph->xincr, ph->array[j][i]);
752                         }
753                 }
754         }
755         hprint(fp, temp);
756         freeharray(temp->array);
757         rfree(temp);
758 }
759 
760 void            hpyprint(FILE * fp, shistogram * ph)
761 {
762         int             i, j;
763         char           *title;
764         shistogram     *temp = NULL;
765 
766         if (ph == NULL)
767         {
768                 return;
769         } else
770         {
771                 title = (char *) ralloc(256 * sizeof(char));
772                 sprintf(title, " Y projection > %s ", ph->title);
773                 temp = hbook1(NHIST - 1, title, ph->ymin, ph->ymax, ph->ybins);
774                 for (j = 0; j < ph->ybins; j++)
775                 {
776                         for (i = 0; i < ph->xbins; i++)
777                         {
778                                 hfill1(temp, ph->ymin + (j + 0.5) * ph->yincr, ph->array[j][i]);
779                         }
780                 }
781         }
782         hprint(fp, temp);
783         freeharray(temp->array);
784         rfree(temp);
785 }
786 
787 double          hdiag(shistogram * ph)
788 {
789         int             j;
790         double          diag = 0.0;
791         if (ph && ph->type == 2)
792         {
793                 for (j = 0; j < ph->ybins; j++)
794                 {
795                         if (j < ph->xbins)
796                         {
797                                 diag += ph->array[j][j];
798                         }
799                 }
800         }
801         return (diag);
802 }
803 
804 void            hintegf(shistogram * ph)
805 {
806         int             i;
807         char           *title;
808 
809         title = (char *) ralloc(256 * sizeof(char));
810 
811         if (ph->type == 1)
812         {
813                 for (i = 1; i < ph->xbins; i++)
814                         ph->array[0][i] += ph->array[0][i - 1];
815                 sprintf(title, " Integrated > %s ", ph->title);
816                 ph->title = title;
817         }
818 }
819 
820 void            hintegb(shistogram * ph)
821 {
822         int             i;
823         char           *title;
824 
825         title = (char *) ralloc(256 * sizeof(char));
826 
827         if (ph->type == 1)
828         {
829                 for (i = ph->xbins; i > 0; i--)
830                         ph->array[0][i - 1] += ph->array[0][i];
831                 sprintf(title, " Integrated < %s ", ph->title);
832                 ph->title = title;
833         }
834 }
835 
836 /* function to differentiate a histogram which has undergone hintegf */
837 void   hdiff(shistogram *ph)
838 {
839   int i;
840 
841   if (ph->type == 1)
842     {
843       for (i = ph->xbins; i > 0; i--)
844         ph->array[0][i] -= ph->array[0][i - 1];
845     }
846 }
847 
848 float hmedian(shistogram *ph)
849 {
850   int i;
851   float median = 0.0, xmin = 0.0;
852   int max = 0.0;
853   float half_max = 0.0;
854   float hist_val1 = 0.0, hist_val2 = 0.0;
855   float fraction = 0.0;
856   int bin_no = 0;
857 
858   if (ph->type == 1)
859     {
860       
861       hintegf(ph);
862       max = ph->contents;
863       xmin = ph->xmin;
864 
865       if (max%2 == 0)
866         {
867           half_max = ((float)max/2)+1;
868         }
869       else
870         half_max = (float)max/2;
871 
872       for (i=0; i<ph->xbins; i++)
873         {
874           hist_val1 = ph->array[0][i];
875           if ((hist_val1 > half_max) && (hist_val2 < half_max))
876             {
877               fraction = (hist_val1 - half_max)/(hist_val1 - hist_val2);
878               bin_no = i;
879             }
880           if (hist_val1 == half_max)
881             {
882               fraction = 0.0;
883               bin_no = i;
884             }
885 
886           hist_val2 = hist_val1;
887         }
888       
889       median = ((bin_no+0.5-fraction)*ph->xincr)+xmin;
890 
891       hdiff(ph);
892     }
893 
894   return(median);
895 }
896 
897 /* essentially the same function as hmedian, but allows the specification of the quartile 
898    instead of imposing the use of 50% */
899 /* NB, quartile should be given in range 0-1.0. ie, 0.5 not 50% */
900 float hquartile(shistogram *ph, float quartile)
901 {
902   int i;
903   float quart_value = 0.0, xmin=0.0;
904   int max = 0.0;
905   float frac_max = 0.0;
906   float hist_val1 = 0.0, hist_val2 = 0.0;
907   float fraction = 0.0;
908   int bin_no = 0;
909 
910   if (ph->type == 1)
911     {
912       
913       hintegf(ph);
914       max = ph->contents;
915       xmin = ph->xmin;
916       if (max%2 == 0)
917         {
918           frac_max = ((float)max*quartile)+1;
919         }
920       else
921         frac_max = (float)max*quartile;
922 
923       for (i=0; i<ph->xbins; i++)
924         {
925           hist_val1 = ph->array[0][i];
926           if ((hist_val1 > frac_max) && (hist_val2 < frac_max))
927             {
928               fraction = (hist_val1 - frac_max)/(hist_val1 - hist_val2);
929               bin_no = i;
930             }
931           if (hist_val1 == frac_max)
932             {
933               fraction = 0.0;
934               bin_no = i;
935             }
936 
937           hist_val2 = hist_val1;
938         }
939       
940       quart_value = ((bin_no+0.5-fraction)*ph->xincr)+xmin;
941 
942       hdiff(ph);
943     }
944 
945   return(quart_value);
946 }
947 
948 
949 void            hsmoof(shistogram * ph)
950 {
951         int             i, j;
952 
953         if (ph->type == 1)
954         {
955                 for (i = 1; i < ph->xbins - 1; i++)
956                         ph->array[0][i - 1] += 0.25 * ph->array[0][i] - 0.125 * ph->array[0][i - 1];
957                 ph->array[0][i + 1] += 0.25 * ph->array[0][i];
958                 ph->array[0][i] = 0.5 * ph->array[0][i];
959         } else if (ph->type == 2)
960         {
961                 for (j = 1; j < ph->ybins - 1; j++)
962                 {
963                         for (i = 1; i < ph->xbins - 1; i++)
964                                 ph->array[j][i] = 0.0625 * ph->array[j - 1][i - 1]
965                                         + 0.0625 * ph->array[j - 1][i]
966                                         + 0.0625 * ph->array[j - 1][i + 1]
967                                         + 0.0625 * ph->array[j][i - 1]
968                                         + 0.5 * ph->array[j][i]
969                                         + 0.0625 * ph->array[j][i + 1]
970                                         + 0.0625 * ph->array[j + 1][i - 1]
971                                         + 0.0625 * ph->array[j + 1][i]
972                                         + 0.0625 * ph->array[j + 1][i + 1];
973                 }
974         }
975 }
976 
977 float           hmax1(shistogram * ph, float *maxx)
978 {
979         int             i;
980         float           hmax = ph->array[0][0];
981         int             imax = 0;
982         double          centre = 0.0;
983         double          h0, h1, hm1;
984 
985         if (ph->type == 1)
986         {
987                 for (i = 0; i < ph->xbins; i++)
988                 {
989                         if (ph->array[0][i] > hmax)
990                         {
991                                 hmax = ph->array[0][i];
992                                 imax = i;
993                         }
994                 }
995         } else
996                 herror("wrong histogram type to hmax1");
997         h0 = hmax;
998         h1 = ph->array[0][imax + 1];
999         hm1 = ph->array[0][imax - 1];
1000         /* assuming parabola fit to 3 centre bins */
1001         centre = (hm1 - h1) / (2 * (h1 + hm1) - 4 * h0);
1002         *maxx = ph->xmin + (imax + 0.5) * ph->xincr;
1003         *maxx = centre * (ph->xincr) + *maxx;
1004 
1005         return (hmax);
1006 }
1007 
1008 float           hnearmax1(shistogram * ph, float *maxx)
1009 {
1010         int             i;
1011         float           hmax = ph->array[0][0];
1012         int             imax = 0;
1013         double          centre = 0.0;
1014         double          h0, h1, hm1, newmaxx, olddiff = FLT_MAX;
1015 
1016         hmax = ph->array[0][0];
1017         imax = 0;
1018 
1019         if (ph->type == 1)
1020         {
1021                 for (i = 1; i < ph->xbins; i++)
1022                 {
1023                         if (ph->array[0][i] > ph->array[0][i - 1] &&
1024                             ph->array[0][i] > ph->array[0][i + 1])
1025                         {
1026                                 if (ph->array[0][i] > hmax)
1027                                 {
1028                                         h0 = ph->array[0][i];
1029                                         h1 = ph->array[0][i + 1];
1030                                         hm1 = ph->array[0][i - 1];
1031                                         /*
1032                                          * assuming parabola fit to 3 centre
1033                                          * bins
1034                                          */
1035                                         centre = (hm1 - h1) / (2 * (h1 + hm1) - 4 * h0);
1036                                         newmaxx = ph->xmin + (i + 0.5) * ph->xincr;
1037                                         newmaxx = centre * (ph->xincr) + newmaxx;
1038                                         if (fabs(newmaxx - *maxx) < olddiff)
1039                                         {
1040                                                 olddiff = fabs(newmaxx - *maxx);
1041                                                 hmax = ph->array[0][i];
1042                                                 imax = i;
1043                                         }
1044                                 }
1045                         }
1046                 }
1047         } else
1048                 herror("wrong histogram type to hmax1");
1049         h0 = hmax;
1050         h1 = ph->array[0][imax + 1];
1051         hm1 = ph->array[0][imax - 1];
1052         /* assuming parabola fit to 3 centre bins */
1053         centre = (hm1 - h1) / (2 * (h1 + hm1) - 4 * h0);
1054         *maxx = ph->xmin + (imax + 0.5) * ph->xincr;
1055         *maxx = centre * (ph->xincr) + *maxx;
1056 
1057         return (hmax);
1058 }
1059 
1060 float           hmax2(shistogram * ph, float *maxx, float *maxy)
1061 {
1062         int             i, j;
1063         float           hmax = ph->array[0][0];
1064         int             imax = 0, jmax = 0;
1065 
1066         if (ph->type == 1)
1067                 herror("wrong histogram type to hmax2");
1068         else
1069         {
1070                 for (j = 0; j < ph->ybins; j++)
1071                 {
1072                         for (i = 0; i < ph->xbins; i++)
1073                         {
1074                                 if (ph->array[j][i] > hmax)
1075                                 {
1076                                         hmax = ph->array[j][i];
1077                                         imax = i;
1078                                         jmax = j;
1079                                 }
1080                         }
1081                 }
1082         }
1083         *maxx = ph->xmin + (imax + 0.5) * ph->xincr;
1084         *maxy = ph->ymin + (jmax + 0.5) * ph->yincr;
1085 
1086         return (hmax);
1087 }
1088 
1089 void            hstore(FILE * fp, shistogram * ph)
1090 /* warning not yet tested */
1091 {
1092         int             i, j;
1093 
1094         fprintf(fp, "%d \n", ph->id);
1095         fprintf(fp, "%s \n", ph->title);
1096         fprintf(fp, "%d \n", ph->type);
1097         fprintf(fp, "%3.2f %3.2f %d \n", ph->xmin, ph->xmax, ph->xbins);
1098         fprintf(fp, "%3.2f %3.2f %d \n", ph->ymin, ph->ymax, ph->ybins);
1099         if (ph->type == 1)
1100         {
1101                 for (i = 0; i < ph->xbins; i++)
1102                         fprintf(fp, "%6.2f %3.2f \n ", ph->xmin + (i + 0.5) * ph->xincr, ph->array[0][i]);
1103         } else if (ph->type == 2)
1104         {
1105                 for (j = 0; j < ph->ybins; j++)
1106                 {
1107                         for (i = 0; i < ph->xbins; i++)
1108                                 fprintf(fp, "%6.2f %6.2f %3.2f \n ",
1109                                         ph->xmin + (i + 0.5) * ph->xincr, ph->ymin + (j + 0.5) * ph->yincr, ph->array[j][i]);
1110                 }
1111                 fprintf(fp, "\n");
1112         }
1113 }
1114 
1115 void            hfetch(FILE * fp)
1116 /* warning not yet tested */
1117 {
1118         shistogram     *ph=NULL;
1119         int             id;
1120         char            title[1024];
1121         int             type;
1122         float           xmin, xmax;
1123         int             xbins;
1124         float           ymin, ymax;
1125         int             ybins;
1126         int             i, j;
1127         float           entry;
1128 
1129 
1130         fscanf(fp, "%d", &id);
1131         fscanf(fp, "%s", title);
1132         fscanf(fp, "%d", &type);
1133         fscanf(fp, "%f %f %d", &xmin, &xmax, &xbins);
1134         fscanf(fp, "%f %f %d", &ymin, &ymax, &ybins);
1135 
1136         if (hist[id] == NULL && type == 1)
1137                 ph = hist[id] = hbook1(id, title, xmin, xmax, xbins);
1138         if (hist[id] == NULL && type == 2)
1139                 ph = hist[id] = hbook2(id, title, xmin, xmax, xbins, ymin, ymax, ybins);
1140         for (j = 0; j < ph->ybins; j++)
1141         {
1142                 for (i = 0; i < ph->xbins; i++)
1143                 {
1144                         if (type == 1)
1145                         {
1146                                 fscanf(fp, "%*f %f", &entry);
1147                                 hfill1(ph, ph->xmin + i * ph->xincr, entry);
1148                         }
1149                         if (type == 2)
1150                         {
1151                                 fscanf(fp, "%*f %*f %f", &entry);
1152                                 hfill2(ph, ph->xmin + i * ph->xincr, ph->ymin + j * ph->yincr, entry);
1153                         }
1154                 }
1155         }
1156 }
1157 
1158 void            hopera()
1159 {
1160 }
1161 
1162 static void     herror(char *str)
1163 {
1164         printf("error **** %s\n", str);
1165         exit(0);
1166 }
1167 
1168 static double   total_chisq(int npar, double *a)
1169 {
1170         double          bh, btot = 0.0;
1171         float           x;
1172         int             j;
1173 
1174         btot = 0.0;
1175         for (x = mph->xmin + 0.5 * (mph->xmax - mph->xmin) / (double) mph->xbins, j = 0;
1176              j < mph->xbins;
1177              x += (mph->xmax - mph->xmin) / (double) mph->xbins, j++)
1178         {
1179                 bh = mph->array[0][j] - mfunc(npar, a, x);
1180                 if (merror_func != NULL)
1181                         bh *= 0.5 * bh / (merror_func(mph, x) * merror_func(mph, x));
1182                 else
1183                         bh *= 0.5 * bh;
1184                 btot += bh;
1185         }
1186         return (btot);
1187 }
1188 
1189 void            hfit(FILE * fp, shistogram * ph, int n, double *a, double *w,
1190                                    double (*fitfunc) (int, double *, float),
1191                                      double (*errfunc) (shistogram *, float))
1192 {
1193         double          min;
1194         double          simplexmin2(int n, double *x, double *lambda,
1195                                    double (*funk) ( /* ??? */ ), void *data,
1196                                    double ftol, void (*text) ( /* ??? */ ));
1197 
1198         mfunc = fitfunc;
1199         merror_func = errfunc;
1200         mph = ph;
1201         min = (float) simplexmin2(n, a, w, total_chisq, NULL, 0.000001, NULL);
1202 }
1203 
1204 double        **herror_analysis(FILE * fp, shistogram * ph, int n, double *a,
1205                                    double (*fitfunc) (int, double *, float),
1206                                    double (*derfunc) (double *, float, int),
1207                                     double (*errfunc) (shistogram *, float),
1208                                                 int exclude)
1209 /* a single data point may be selected for exclusion */
1210 {
1211         double        **c=NULL, **cinv=NULL, *placebo=NULL;
1212         int             i, j, k;
1213         double          d, x;
1214 
1215         c = maketarray(n, n, 0.0);
1216         for (i = 0; i < n; i++)
1217         {
1218                 for (j = 0; j < n; j++)
1219                 {
1220                         c[i][j] = 0.0;
1221                         for (x = ph->xmin + ph->xincr / 2.0, k = 0; x < ph->xmax; x += ph->xincr, k++)
1222                         {
1223                                 if (k != exclude)
1224                                 {
1225                                         d = (double) (*errfunc) (ph, x);
1226                                         d *= d;
1227                                         c[i][j] += (double) (*derfunc) (a, x, i) * (double) (*derfunc) (a, x, j) / d;
1228                                 }
1229                         }
1230                 }
1231         }
1232 
1233         cinv = hmatinv(n, a, c, placebo);
1234         freetarray(c);
1235 
1236         if (fp != (FILE *) 0)
1237         {
1238                 for (i = 0; i < n; i++)
1239                 {
1240                         fprintf(fp, "     parameter %d = %f +- %f\n", i, a[i], sqrt(cinv[i][i]));
1241                 }
1242         }
1243         return (cinv);
1244 }
1245 
1246 void            hsuper(shistogram * ph, int n, double (*func) (int, double *, float), double *a)
1247 {
1248         int             i;
1249 
1250         if (ph->npar != n && ph->par != NULL)
1251         {
1252                 rfree(ph->par);
1253                 ph->par = NULL;
1254         }
1255         if (ph->par == NULL)
1256                 ph->par = (double *) ralloc(n * sizeof(double));
1257         ph->shFuncSuper = func;
1258         ph->npar = n;
1259         for (i = 0; i < n; i++)
1260                 ph->par[i] = a[i];
1261 }
1262 
1263 void            hfit_gauss(FILE * fp, shistogram * ph)
1264 {
1265         double          a[3];
1266         double          w[3];
1267         int             i;
1268 
1269         if (ph == NULL || fp == NULL)
1270                 return;
1271 
1272         a[1] = ph->mean;
1273         a[2] = (double) sqrt((double) (ph->mean2 - ph->mean * ph->mean));
1274 
1275         a[0] = 0.0;
1276 
1277         for (i = 0; i < ph->xbins; i++)
1278                 if (ph->array[0][i] > a[0])
1279                         a[0] = ph->array[0][i];
1280 
1281         w[0] = a[0];
1282         w[1] = a[2];
1283         w[2] = a[2];
1284 
1285         if (a[2] == 0.0)
1286                 return;
1287 
1288         hfit(fp, ph, 3, a, w, hgaussian, NULL);
1289 }
1290 
1291 double          herrdum(shistogram * ph, float x)
1292 {
1293         double          d;
1294         double          y;
1295 
1296         y = hfill1(ph, x, 0.0);
1297 
1298         if (y <= 0.0)
1299                 y = 1.0;
1300         d = sqrt(y);
1301 
1302         return (d);
1303 }
1304 
1305 double          hgaussian(int n, double *a, float x)
1306 {
1307         return (a[0] * exp(-(a[1] - x) * (a[1] - x) / (2.0 * a[2] * a[2])));
1308 }
1309 
1310 double          hgaussder(double *a, float x, int i)
1311 {
1312         double          y;
1313         y = hgaussian(3, a, x);
1314         if (i == 0)
1315                 return (exp(-(a[1] - x) * (a[1] - x) / (2.0 * a[2] * a[2])));
1316 
1317         if (i == 1)
1318                 return ((x - a[1]) * y / (a[2] * a[2]));
1319 
1320         if (i == 2)
1321                 return ((x - a[1]) * (x - a[1]) / (a[2] * a[2] * a[2]) * y);
1322 
1323         return (0.0);           /* should never be executed */
1324 }
1325 
1326 double        **hmatinv(int n, double *params, double **a, double *nresult)
1327 {
1328         double        **y, **p, **q, *d, *nparams;
1329         double          temp;
1330         /*
1331             double condition = 0.000001;
1332         */
1333         double          condition = 1.0;
1334         int             i, j, k;
1335 
1336         y = maketarray(n, n, 0.0);
1337         p = maketarray(n, n, 0.0);
1338         q = maketarray(n, n, 0.0);
1339         d = (double *) ralloc(n * sizeof(double));
1340         nparams = (double *) ralloc(n * sizeof(double));
1341         svd(n, n, a, p, d, q);
1342 
1343         /*
1344             wmax = 0.0;
1345             for (j = 0; j < n; j++)
1346                 if (d[j] > wmax)
1347                     wmax = d[j];
1348             thresh = condition * wmax;
1349             for (j = 0; j < n; j++)
1350                 if (d[j] < thresh)
1351                     d[j] / thresh*thresh;
1352                 else
1353                     d[j] = 1.0/d[j];
1354         */
1355         /* compute rotated parameter set */
1356         for (i = 0; i < n; i++)
1357         {
1358                 nparams[i] = 0.0;
1359                 for (j = 0; j < n; j++)
1360                 {
1361                         nparams[i] += params[j] * p[j][i];
1362                 }
1363                 nparams[i] *= nparams[i];
1364         }
1365         /* test new parameter set for significant non-zero */
1366         *nresult = 0.0;
1367         for (j = 0; j < n; j++)
1368         {
1369                 temp = d[j];
1370                 if (d[j] * nparams[j] < condition)
1371                         d[j] *= nparams[j] * nparams[j] / (condition * condition);
1372                 else
1373                         d[j] = 1.0 / d[j];
1374                 *nresult += temp * d[j];
1375         }
1376 
1377         for (i = 0; i < n; i++)
1378         {
1379                 for (j = 0; j < n; j++)
1380                 {
1381                         for (k = 0; k < n; k++)
1382                         {
1383                                 y[j][k] += d[i] * p[j][i] * q[k][i];
1384                         }
1385                 }
1386         }
1387 
1388         freetarray(p);
1389         freetarray(q);
1390         rfree(d);
1391         return (y);
1392 }
1393 
1394 
1395 /*
1396         a.lacey additions
1397 */
1398 
1399 
1400 
1401 void            hbin_dump_plain(FILE * out, shistogram * ph)
1402 {
1403         int             i;
1404         float           c;
1405 
1406         if ((ph == NULL) || (out == NULL))
1407                 return;
1408 
1409         c = ph->xmin;
1410         for (i = 0; i < ph->xbins; i++)
1411         {
1412                 fprintf(out, "%f\t%f\n", c, ph->array[0][i]);
1413                 c += ph->xincr;
1414         }
1415 }
1416 

~ [ source navigation ] ~ [ diff markup ] ~ [ identifier search ] ~ [ freetext search ] ~ [ file search ] ~

This page was automatically generated by the LXR engine.
Visit the LXR main site for more information.