/* dst_to_src... c                             Steve Mann, 1993 */
/* mapping for use by dechirp2.c, rechirp2.c, cement.f, ... */

/* global to this file */
static int in_vdim, in_hdim;
static int out_vdim, out_hdim;
static double a_inv,b_inv,c_inv,d_inv,e_inv,f_inv,g_inv,h_inv;

/* converts m1..n4 to p-chirplet parameters; returns nothing */
/* sets global variables for this file */
void dst_to_src_once_per_image(M_in,N_in,
                               M_out,N_out,
                               c1,c2,c3,c4,c5,c6,c7,c8) /* chirp parameters */
int M_in, N_in;
int M_out, N_out;
double c1,c2,c3,c4,c5,c6,c7,c8;
{
      in_vdim=M_in;
      in_hdim=N_in;
      out_vdim=M_out;
      out_hdim=N_out;

      /* pCHIRP CODE */

      a_inv = c1;  /* replaced 1500 lines or so with this identity map */
      b_inv = c2;
      c_inv = c3;
      d_inv = c4;
      e_inv = c5;
      f_inv = c6;
      g_inv = c7;
      h_inv = c8;
}

/* takes destination pixel coordinates [m,n] and
   returns source sub-pixel coordinates [m_out, n_out] */
void dst_to_src_once_per_pixel(m,n,m_in,n_in)
int m, n;
double *m_in, *n_in;
{
  double x_in, y_in;
  double x_d, y_d;
          /* normalize input scan, x_in and y_in, to [0,1) */
          x_in=(double)m/(double)out_vdim;
          y_in=(double)n/(double)out_hdim;

          /* if( DEBUG ) fprintf(stderr,"x_in=%f  y_in=%f\n",x_in,y_in); */
          /* incomplete Lorentz transformations; replaced i with 1.0 */
          /* mapping from input vars (x_in,y_in) to input vars (x,y) */
          x_d = (a_inv*x_in + b_inv*y_in + c_inv)/
                (g_inv*x_in + h_inv*y_in + 1.0);
          y_d = (d_inv*x_in + e_inv*y_in + f_inv)/
                (g_inv*x_in + h_inv*y_in + 1.0);
          *m_in = (x_d*(double)in_vdim);
          *n_in = (y_d*(double)in_hdim);
}

void dst_to_src_once_per_pixel_noscaling(m,n,m_in,n_in)
int m, n;
double *m_in, *n_in;
{
  double x_in, y_in;
  double x_d, y_d;
          /* normalize input scan, x_in and y_in, to [0,1) */
          x_in=(double)m/(double)in_vdim;
          y_in=(double)n/(double)in_hdim;

          /* if( DEBUG ) fprintf(stderr,"x_in=%f  y_in=%f\n",x_in,y_in); */
          /* incomplete Lorentz transformations; replaced i with 1.0 */
          /* mapping from input vars (x_in,y_in) to input vars (x,y) */
          x_d = (a_inv*x_in + b_inv*y_in + c_inv)/
                (g_inv*x_in + h_inv*y_in + 1.0);
          y_d = (d_inv*x_in + e_inv*y_in + f_inv)/
                (g_inv*x_in + h_inv*y_in + 1.0);
          *m_in = (x_d*(double)in_vdim);
          *n_in = (y_d*(double)in_hdim);
}