/*
   File:        contours.c
   Author:      Andrew W. Moore
   Created:     Fri May 14 16:32:22 EDT 1993
   Description: Drawing contours of functions and 2d arrays

   Copyright (C) 1992, Andrew W. Moore
*/

#include <stdio.h>
#include <math.h>
#include "ambs.h"      /* Very basic operations */
#include "amgr.h"      /* Basic (0,512)x(0,512) Graphics window */
#include "amma.h"      /* Fast, non-fragmenting, memory management */
#include "amar.h"      /* Obvious operations on 1-d arrays */
#include "geom.h"      /* Simple 2-d geometry structures and operations */
#include "contours.h"  /* Drawing contours of functions and 2d arrays */

condata *make_condata(num_x_parts,num_y_parts)
int num_x_parts,num_y_parts;
{
  condata *new = AM_MALLOC(condata);
  int j;

  new -> num_x_parts = num_x_parts;
  new -> num_y_parts = num_y_parts;
  new -> data = AM_MALLOC_ARRAY(float_ptr,num_y_parts+1);

  for ( j = 0 ; j <= num_y_parts ; j++ )
  {
    new->data[j] = AM_MALLOC_ARRAY(float,num_x_parts+1);
    set_floats_constant(new->data[j],num_x_parts+1,0.0);
  }

  return(new);
}

void free_condata(cd)
condata *cd;
{
  int j;
  for ( j = 0 ; j < cd -> num_y_parts ; j++ )
    am_free((char *) (cd->data[j]),(cd->num_x_parts+1) * sizeof(float));

  am_free((char *) cd->data,(cd->num_y_parts+1) * sizeof(float));
  am_free((char *) cd,sizeof(condata));
}

float condata_ref(cd,i,j)
condata *cd;
int i,j;
{
  if ( i < 0 || i > cd->num_x_parts || j < 0 || j > cd->num_y_parts )
    my_error("otyihjyhtoitoxgac");
  else
  {
    float *fp = cd->data[j];
    return(fp[i]);
  }
}

void condata_update(cd,i,j,v)
condata *cd;
int i,j;
float v;
{
  if ( Verbosity > 100.0 )
    printf("condata_update(cd,%d,%d,%g)\n",i,j,v);

  if ( i < 0 || i > cd->num_x_parts || j < 0 || j > cd->num_y_parts )
    my_error("otyihjyhtoitoxgac");
  else
  {
    float *fp = cd->data[j];
    fp[i] = v;
  }
}

void find_lambda_intersection(ct,n1,n2,height,p)
contrig *ct;
int n1,n2;
float height;
point *p;
/*
   PRE: the value "height" lies between
        the values at corners n1 and n2 of the triangle.

   POST: Performing linear interpolation betwen them, p contains the 2-d
         point on the triangle edge which actually should equal "height"

   ALGORITHM:
     p = c[n1] + lambda (c[n2] - c[n1])
     height = v[n1] + lambda(v[n2] - v[n1])
       ...know everything but lambda and p. Solve for lambda then p.
     
*/
{
  float lambda = (height - ct->values[n1]) / (ct->values[n2] - ct->values[n1]);
  subtract_points(&(ct->corners[n2]),&(ct->corners[n1]),p);
  p->x *= lambda;
  p->y *= lambda;
  add_points(&(ct->corners[n1]),p,p);
}

lines *find_contrig(ct,height,ls)
contrig *ct;
float height;
lines *ls;
{
  bool under0 = height < ct->values[0];
  bool under1 = height < ct->values[1];
  bool under2 = height < ct->values[2];

  if ( under0 != under1 || under1 != under2 )
  {
    int odd;
    int same_a,same_b;
    line l;

    if ( under0 == under1 )
      odd = 2;
    else if ( under1 == under2 )
      odd = 0;
    else
      odd = 1;

    same_a = (odd + 1) % THREE;
    same_b = (odd + 2) % THREE;

    find_lambda_intersection(ct,odd,same_a,height,&l.start);    
    find_lambda_intersection(ct,odd,same_b,height,&l.end);

    ls = add_to_lines(&l,ls);
  }
  return(ls);
}

float caxis_value(cax,i)
caxis *cax;
int i;
{
  return(cax->lo + i * (cax->hi - cax->lo)/cax -> parts);
}

lines *find_contrig_contours(cen,ct,ls)
conenv *cen;
contrig *ct;
lines *ls;
{
  int k;

  if ( Verbosity > 100.0 )
  {
    int i;
    printf("contrig contours: ");
    for ( i = 0 ; i < THREE ; i++ )
   printf(" [(%g,%g) -> %g] ",ct->corners[i].x,ct->corners[i].y,ct->values[i]);
    printf("\n");

    wait_for_key();
  }

  for ( k = 0 ; k <= cen->conaxes->z.parts ; k++ )
    ls = find_contrig(ct,caxis_value(&cen->conaxes->z,k),ls);
  return(ls);
}

void contrig_corner_from_data(cen,i,j,n,ct)
conenv *cen;
int i,j,n;
contrig *ct;
{
  ct->corners[n].x = caxis_value(&cen->conaxes->x,i);
  ct->corners[n].y = caxis_value(&cen->conaxes->y,j);
  ct->values[n] = condata_ref(cen->condata,i,j);
}

lines *find_partition_contours(cen,i,j,ls)
conenv *cen;
int i,j;
lines *ls;
{
  contrig ct;
  contrig_corner_from_data(cen, i   , j   , 0 ,&ct);
  contrig_corner_from_data(cen, i+1 , j+1 , 1 ,&ct);

  contrig_corner_from_data(cen, i+1 , j   , 2 ,&ct);
  ls = find_contrig_contours(cen,&ct,ls);

  contrig_corner_from_data(cen, i   , j+1 , 2 ,&ct);
  ls = find_contrig_contours(cen,&ct,ls);
  return(ls);
}

lines *find_conenv_contours(cen,ls)
conenv *cen;
lines *ls;
/*
   PRE: cen->conaxes and cen->condata defined and in agreement. 
*/
{
  int i,j;

  for ( i = 0 ; i < cen->conaxes->x.parts ; i++ )
    for ( j = 0 ; j < cen->conaxes->y.parts ; j++ )
      ls = find_partition_contours(cen,i,j,ls);

  return(ls);
}

void init_caxis(ax,lo,hi,parts)
caxis *ax;
float lo,hi;
int parts;
{
  ax->lo = lo;
  ax->hi = hi;
  ax->parts = parts;
}

condata *make_condata_from_function(cax,height_fn,height_fn_data)
conaxes *cax;
float (*height_fn)();
char *height_fn_data;
/*
    height_fn : (float) x , (float) y , (char *) data -----> (float) height
*/
{
  condata *cd = make_condata(cax->x.parts,cax->y.parts);
  int i,j;
  for ( i = 0 ; i <= cax->x.parts ; i++ )
    for ( j = 0 ; j <= cax->y.parts ; j++ )
    {
      float x = caxis_value(&cax->x,i);
      float y = caxis_value(&cax->y,j);
      float height = (*height_fn)(x,y,height_fn_data);
      condata_update(cd,i,j,height);

      if ( Verbosity > 200.0 )
      {
        int ti,tj;
        for ( ti = 0 ; ti <= i ; ti++ )
         for ( tj = 0 ; tj <= cd->num_y_parts && (tj<=j || ti<i) ; tj++ )
          printf("  condata_ref(cd,%d,%d) = %g\n",ti,tj,condata_ref(cd,ti,tj));
        wait_for_key();
      }
    }

  return(cd);
}

lines *find_scalefree_contours(cd,z_caxis)
condata *cd;
caxis *z_caxis;
{
  conaxes cax;
  conenv cen;
  lines *ls;

  init_caxis(&cax.x,0.0,512.0,cd->num_x_parts);
  init_caxis(&cax.y,0.0,512.0,cd->num_y_parts);
  init_caxis(&cax.z,z_caxis->lo,z_caxis->hi,z_caxis->parts);

  cen.conaxes = &cax;
  cen.condata = cd;
  
  ls = find_conenv_contours(&cen,(lines *)NULL);
  return(ls);
}


lines *find_function_contours(cax,height_fn,height_fn_data)
conaxes *cax;
float (*height_fn)();
char *height_fn_data;
/*
    height_fn : (float) x , (float) y , (char *) data -----> (float) height
*/
{
  condata *cd = make_condata_from_function(cax,height_fn,height_fn_data);
  lines *ls = find_scalefree_contours(cd,&cax->z);
  free_condata(cd);
  return(ls);
}

void draw_data_contours(cd,z_caxis)
condata *cd;
caxis *z_caxis;
{
  lines *ls = find_scalefree_contours(cd,z_caxis);
  draw_lines(ls);
  free_lines(ls);
}

void draw_function_contours(cax,height_fn,height_fn_data)
conaxes *cax;
float (*height_fn)();
char *height_fn_data;
/*
    height_fn : (float) x , (float) y , (char *) data -----> (float) height
*/
{
  condata *cd = make_condata_from_function(cax,height_fn,height_fn_data);
  draw_data_contours(cd,&cax->z);
  free_condata(cd);
}