#include <math.h>
 /*
 /* Math.h supplies double precision exp(), pow(), and sinh() functions.
  * If your math.h doesn't supply sinh(), use sinh(x)=(exp(x)+exp(-x))/2.
  */
double pcal_div(double n, double d)
{
   return floor((double)n/(double)d);
}

/*
 * pcal_make_lut
 *
 * Returns a lookup table for recovering physical values
 * from pCAL-encoded PNG samples
 *
 * Glenn Randers-Pehrson, US Army Research Laboratory, January 1997
 * glennrp@arl.mil or randeg@alumni.rpi.edu
 *
 * This code still may need some bullet-proofing in case
 * intermediate terms of the products overflow.
 *
 */
int pcal_make_lut (float *physical_value, unsigned int m, int equation_type,
                   long x0, long x1, float *p)

/* returns 0 (success)
 *        -1 (error, x0==x1)
 *        -2 (unknown equation type)
 * input:
 *         m: PNG sample depth
 *         equation_type: from pCAL chunk
 *         x0, x1: stored sample to original sample mapping from pCAL chunk
 *         p[]: equation parameters, from pCAL chunk
 * output:
 *         physical_value[0..m] (caller must allocate space for the array)
 */

{
    long sample, osample;
    float fd;
    double d, dm, two;

    two=2.;
    d=x1-x0;
    dm=m;
    fd=d;

    if(x0 == 0 && x1 == m) {

    /* The usual case, where stored_samples are original
     * samples. It is ok to use regular integer division
     * because the operands are never negative.
     */
        if(equation_type == 0){
            for (sample=0; sample<=m; sample++) physical_value[sample] =
                   p[0] + p[1]*sample/fd;
        }
        else if(equation_type == 1){
            for (sample=0; sample<=m; sample++) physical_value[sample] =
                   p[0] + p[1]*exp(p[2]*(sample/fd));
        }
        else if(equation_type == 2){
            for (sample=0; sample<=m; sample++) physical_value[sample] =
                   p[0] + p[1]*pow(p[2],(float)(sample/fd));
        }
        else if(equation_type == 3){
            for (sample=0; sample<=m; sample++) physical_value[sample] =
                  p[0] + p[1]*sinh(p[2]*(sample-p[3])/fd);
        }
        else return (-2); /* ERROR, unknown equation type */
    }
    else if (x1 != x0) {
    /*
     * It is not ok to use regular integer division because
     * the operands might be negative; also original_samples
     * must be calculated.
     */

     /*
      *   original_sample = 
      *      (stored_sample * (X1 - X0) + M/2) / M + X0
      */

        if(equation_type == 0){
            for (sample=0; sample<=m; sample++){
               osample=pcal_div(sample*d+pcal_div(dm,two),dm) + x0;
               physical_value[sample] = p[0] + p[1]*osample/fd;
            }
        }
        else if(equation_type == 1){
            for (sample=0; sample<=m; sample++){
               osample=pcal_div(sample*d+pcal_div(dm,two),dm) + x0;
               physical_value[sample] = p[0] + p[1]*exp(p[2]*osample/fd);
            }
        }
        else if(equation_type == 2){
            for (sample=0; sample<=m; sample++){
               osample=pcal_div(sample*d+pcal_div(dm,two),dm) + x0;
               physical_value[sample] = p[0] + p[1]*pow(p[2],osample/fd);
            }
        }
        else if(equation_type == 3){
            for (sample=0; sample<=m; sample++){
               osample=pcal_div(sample*d+pcal_div(dm,two),dm) + x0;
               physical_value[sample] =
                  p[0] + p[1]*sinh(p[2]*(osample-p[3])/fd);
            }
        }
        else return (-2); /* ERROR, unknown equation type */
    }
    else return (-1); /* ERROR, X0 == X1 */

    return (0);
}
unsigned short limit(long low, float x, long high)
{
    if (low < high){
       if(x < low)return low;
       else if (x > high) return high;
       else return x;
    }
    else {
       if(x < high)return high;
       else if( x > low) return low;
       else return x;
    }
}
/*
 * pcal_encode
 *
 * encodes an array of physical values as pCAL-encoded PNG samples
 *
 * Glenn Randers-Pehrson, US Army Research Laboratory, January 1997
 * glennrp@arl.mil or randeg@alumni.rpi.edu
 *
 */
#include <math.h>
 /* Math.h supplies double precision exp(), pow(), and sinh() functions.
  * If your math.h doesn't supply sinh(), use sinh(x)=(exp(x)+exp(-x))/2.
  */
int pcal_encode (unsigned short *stored_sample, long n, float *physical_value, 
 unsigned int m, int equation_type, long x0, long x1, float *p)

/* returns 0 (success)
 *        -1 (error, x0==x1)
 *        -2 (unknown equation type)
 * input:
 *         n: number of samples
 *         physical_value[0..n-1]
 *         m: PNG sample depth
 *         equation_type: from pCAL chunk
 *         x0, x1: stored sample to original sample mapping from pCAL chunk
 *         p[]: equation parameters, from pCAL chunk
 * output:
 *         stored_samples[0..n-1] (caller must allocate space for the array)
 */

{
    long isample, osample;
    double d,dm;
    float fd;
    long k;

    d=x1-x0;
    fd=d;
    dm=m;

    if(x0 == 0 && x1 == m) {

    /* The usual case, where stored_samples are original
     * samples. It is ok to use regular integer division
     * because the operands are never negative.
     */
        if(equation_type == 0){
            for (k=0; k<n; k++) stored_sample[k]=
              limit(0, .5+(d*(physical_value[k] - p[0])/p[1]), m);
        }
        else if(equation_type == 1){
            for (k=0; k<n; k++) stored_sample[k]=
              limit(0, .5+d*(log(physical_value[k] - p[0])/p[1])/p[2], m);
        }
        else if(equation_type == 2){
            float factor;
            factor=fd/log(p[2]);
            for (k=0; k<n; k++) stored_sample[k]=
              limit(0, .5+log((physical_value[k]-p[0])/p[1])*factor, m);
        }
        else if(equation_type == 3){
            for (k=0; k<n; k++) stored_sample[k]=
               stored_sample[k]=limit(0,
                    .5+p[3]+fd*asinh((physical_value[k]-p[0])/p[1])/p[2], m);
        }
        else return (-2); /* ERROR, unknown equation type */
    }
    else if (x1 != x0) {
        long two;
        two=2.0;
    /*
     * It is not ok to use regular integer division because
     * the operands might be negative; also original_samples
     * must be calculated.
     *  ((original_sample - X0) * M + (X1 - X0) / 2) / (X1 - X0)
     */
        if(equation_type == 0){
            for (k=0; k<n; k++) {
               isample=.5+d*(physical_value[k] - p[0])/p[1];
               osample=limit(x0, isample, x1);
               stored_sample[k]= pcal_div((osample-x0)*dm+pcal_div(d,two),d);
            }
        }
        else if(equation_type == 1){
            for (k=0; k<n; k++) {
               isample= .5+d*(log(physical_value[k] - p[0])/p[1])/p[2];
               osample=limit(x0, isample, x1);
               stored_sample[k]= pcal_div((osample-x0)*dm+pcal_div(d,two),d);
            }
        }
        else if(equation_type == 2){
            float factor;
            factor=fd/log(p[2]);
            for (k=0; k<n; k++) {
               isample=.5+log((physical_value[k]-p[0])/p[1])*factor;
               osample=limit(x0, isample, x1);
               stored_sample[k]= pcal_div((osample-x0)*dm+pcal_div(d,two),d);
            }
        }
        else if(equation_type == 3){
            for (k=0; k<n; k++) {
               isample= .5+p[3]+fd*asinh((physical_value[k]-p[0])/p[1])/p[2];
               osample=limit(x0, isample, x1);
               stored_sample[k]= pcal_div((osample-x0)*dm+pcal_div(d,two),d);
            }
        }
        else return (-2); /* ERROR, unknown equation type */
    }
    else return (-1); /* ERROR, x0 == x1 */

    return (0);
}
int pcal_print( float *physical_value, unsigned int m)
{
    long sample;
    unsigned int skip;
    skip=m/8+1;
    if(m < 100){
        for (sample=0;sample<=m;sample+=skip)
            printf("  p[%2d] = %8e\n",sample,physical_value[sample]);
        if(sample < m+skip)
            printf("  p[%2d] = %8e\n",m,physical_value[m]);
        }
    else
        {
        for (sample=0;sample<=m;sample+=skip)
            printf("  p[%5d] = %8e\n",sample,physical_value[sample]);
        if(sample < m+skip)
            printf("  p[%5d] = %8e\n",m,physical_value[m]);
        printf("\n");
        if(m = 65535)
            for (sample=32765;sample<=32770;sample++)
                printf("  p[%5d] = %8e\n",sample,physical_value[sample]);
    }
    return (0);
}
#include <math.h>
main()
{
    float pval[65536], p[4];
    float original_pval[64];
    unsigned short samples[64];
    long x0,x1;
    unsigned int m;
    int equation_type;
    int i;
    long n=32;

    for (i=0; i<n; i+=4){
          float exponent;
          exponent=(i-4)/4.0;
          original_pval[i]=pow(10.,exponent);
          original_pval[i+1]=2.*original_pval[i];
          original_pval[i+2]=5.*original_pval[i];
          original_pval[i+3]=9.*original_pval[i];
    }
     
    m=65535;
    x0=0;
    x1=64000;

    equation_type=0;
    p[0]=0;
    p[1]=1000000.;

    printf("\n test 0: 0-1,000,000, linear\n");
    pcal_encode(samples,n,original_pval,m,equation_type,x0,x1,p);
    pcal_make_lut(pval,m,equation_type,x0,x1,p);
    pcal_print(pval, 7);
    pcal_show (32,original_pval,samples,pval);
    
    equation_type=1;
    p[1]=1.;
    p[2]=log(1000000.);

    printf("\n test 1: 1-1,000,000, exp()\n");
    pcal_encode(samples,n,original_pval,m,equation_type,x0,x1,p);
    pcal_make_lut(pval,m,equation_type,x0,x1,p);
    pcal_print(pval, 7);
    pcal_show (32,original_pval,samples,pval);

    equation_type=2;
    p[1]=1.;
    p[2]=1000000.;

    printf("\n test 2: 1-1,000,000, pow()\n");
    pcal_encode(samples,n,original_pval,m,equation_type,x0,x1,p);
    pcal_make_lut(pval,m,equation_type,x0,x1,p);
    pcal_print(pval, 7);
    pcal_show (32,original_pval,samples,pval);
    
    equation_type=3;
    p[0] = 0.0;
    p[1] = 1.0e-29;
    p[2] = 280.0;
    p[3] = 32767.0;

    printf("\n test 3: full range of 32-bit floats, sinh()\n");
    pcal_encode(samples,n,original_pval,m,equation_type,x0,x1,p);
    pcal_make_lut(pval,m,equation_type,x0,x1,p);
    pcal_print(pval, 7);
    pcal_show (32,original_pval,samples,pval);
}
pcal_show(long n,float *original_pval, unsigned short *samples,
          float *physical_value)
{
    float decoded_value, error, error_percent;
    long i;

    printf(
"\n    original    stored      decoded       total     percent \n");
    printf(
  "    physical    sample        value       error       error \n");
    printf(
  "   ---------    ------     --------     -------     ------- \n");
    for (i=0; i<n; i++){
        long percent;
        decoded_value=physical_value[samples[i]];
        error=decoded_value - original_pval[i];
        percent=1000000.*error/original_pval[i];
        error_percent=percent/10000.;
        printf(" %11.4g %10d %11.4g %11.4g %11.4g\n",
              original_pval[i],samples[i],decoded_value,
              error, error_percent);
     }
}