#include <float_image.h>
#include <pst_imgsys.h>
#include <pst_slope_map.h>
#include <pst_weight_map.h>
#include <pst_height_map.h>

void pst_slope_and_weight_map_shrink(float_image_t *IG, float_image_t *IW, float_image_t** SG, float_image_t** SW);
void weighted_average(double w0, double v0, double w1, double v1,double vardelta, double *v, double *w);
/* Given two estimates {v0 = vm - delta/2 + e0} and {v1 = vm + delta/2 + e1} of a quantity {vm},
  contaminated by a common error {delta} and indepedent errors {e0,e1}, returns in {*v} the most reliable
  estimate of {vm}.  Assumes that {delta,e0,e1} are Gaussian variables with variances {vardelta, 1/w0, and 1/w1},
  respectively. It works even when these variances are 0 or {+INF}. */

imgsys_t* pst_slope_map_build_integration_system(float_image_t *G, float_image_t *GW, bool_t full, float_image_t *ZW);

void pst_slope_map_to_depth_map_recursive
  ( float_image_t *IG, 
    float_image_t *IW, 
    int level,
    bool_t newStyle,
    long int maxIter,
    long int minMapSize,
    double convTol,
    float_image_t **OZ,
    float_image_t **OW,
    bool_t verbose,
    long int reportIter,
    pst_slope_map_report_proc_t *reportData,
    pst_imgsys_report_proc_t *reportSys,
    pst_height_map_report_proc_t *reportHeights
  );
  

imgsys_t *pst_imgsys_from_eqs(int NX, int NY, long int N, imgsys_equation_t *eq, long int *ix, int *col, int *row);
  /* Creates a new linear system with {N} equations on {N} unknowns
    for an image with {NX} columns and {NY} rows, using the given
    index mapping tables {ix,col,row} and the given equation vector {eq}. Does not allocate
    any tables and does not initialize or modify the equations. */

void pst_compute_edge_slope_and_weight
  ( float_image_t *G,
    float_image_t *GW,
    int axis,
    int x,int y,
    double *d,
    double *w
  );
  /* Computes the derivative {*d} and weight{*w} along the edge with direction
     {axis} (0 = X , 1 = Y) whose origin is the vertex {(x,y)}. */
  
void pst_slope_map_solve_system
  ( imgsys_t *S, 
    float_image_t *OZ, 
    float_image_t *OW,
    long int maxIter, 
    double convTol, 
    int para, 
    int szero, 
    bool_t verbose, 
    int level, 
    long int reportIter, 
    pst_height_map_report_proc_t *reportHeights
  );
  /* 
    Solves system {S} by Gauss-Seidel iteration, and stores the solution
    into the height map {OZ}.  On input {OZ} should contain the initial
    guess.
      
    The image {OZ} must be single-channel. The procedure assumes that the
    unknowns of {S} are related to the samples of {OZ}, as defined by
    the tables {ix[0..NX*NY-1]}, {col[0..S->N-1]}, and {row[0..S->N-1]}.
    
    The image {OW} is supposed to be a weight map with the same dimensions
    as {OZ}. It is merely passed to {reportHeights}.

    If {reportHeights} is not null and {reportIter} is not zero, will call
    {reportZ(level,iter,change,final,OZ,OW)} before the first iteration, after
    the last iteration, and after any number of iterations that is
    divisible by {reportIter}. */

void pst_interpolate_two_samples
  (  float_image_t* I, float_image_t* W,
     int c,
     int x0, int y0,
     int x1, int y1,
     double *v, double* w
   );
  
  
void pst_interpolate_four_samples
  (  float_image_t* I, float_image_t* W,
     int c,
     int xm, int ym,
     int x0, int y0,
     int x1, int y1,
     int xp, int yp,
     double *v, double* w
   );
  
void pst_interpolate_four_points
  (  double vm, double wm,
     double v0, double w0,
     double v1, double w1,
     double vp, double wp,
     double *v, double *w
  );	
  
long int* pst_imgsys_sort_equations(imgsys_t *S, float_image_t *OW);