#include "VsipTypes.h"

/****************************************************************************
 *
 *  Function    : Vsip8x8InvDiscreteCosineTransform
 *
 *  Description : Compute the 8 x 8 Inverse Discrete Cosine Transform of an
 *                Image. The Image must have size multiple of 8x8
 *
 *  Parameters  : src - input image object
 *
 *  Note        : The IntType is computed using the IOtype and assuming
 *                the actual alpha values (N = 8)
 *
 ****************************************************************************/
$generic

   $Vsip8x8InvDiscreteCosineTransform = 'Vsip8x8InvDiscreteCosineTransformByte',
        $InType = 'int32', $OutType = 'uint8', $InterType = 'fix32.16';
   $Vsip8x8InvDiscreteCosineTransform = 'Vsip8x8InvDiscreteCosineTransformInt',
        $InType = 'int32', $OutType = 'int32', $InterType = 'fix32.16';

$in

$OutType [:,:] $Vsip8x8InvDiscreteCosineTransform($InType src[:,:])
{

  // ***** recover the src size *****
  VsipIndexesType rs, VsipIndexesType cs = extents(src);

 // ***** Test size of the input image (src) (see header) *****
 assert(((rs > 0) && (cs > 0)),
         "ERROR: Image can not have zero dimention. (",rs,"x",cs,")\n");
 assert((rs % 8 == 0) && (cs % 8 == 0),
        "ERROR: Image must be multiple of 8x8.(",rs,"x",cs,")\n");

 // ***** set alpha. alpha(0) = sqrt(1/N).          *****
 // ***** alpha(u) = sqrt(2/N). In this case N = 8  *****
 fix11.10 alpha[8]= {0.35355339, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5};

 fix11.10 alpha2[:,:] =
      for elem1 in alpha cross elem2 in alpha
      return (array(elem1*elem2));

 // ***** Table of precomputed values used inside the convolution      *****
 // ***** cos((2 * i + 1) * u * PI / (2 * N)) for i, u = 0, ..., 7     *****
 // ***** Computed using PI = 3.14159265 and N = 8                     *****
 fix12.10 costab[8,8] = { {  1.00000000,  0.98078528,  0.92387953,  0.83146961,
                             0.70710678,  0.55557023,  0.38268343,  0.19509032},
                          {  1.00000000,  0.83146961,  0.38268343, -0.19509032,
                            -0.70710678, -0.98078528, -0.92387953, -0.55557024},
                          {  1.00000000,  0.55557023, -0.38268343, -0.98078528,
                            -0.70710678,  0.19509032,  0.92387953,  0.83146962},
                          {  1.00000000,  0.19509032, -0.92387953, -0.55557024,
                             0.70710678,  0.83146962, -0.38268342, -0.98078528},
                          {  1.00000000, -0.19509032, -0.92387953,  0.55557023,
                             0.70710679, -0.83146961, -0.38268344,  0.98078528},
                          {  1.00000000, -0.55557023, -0.38268344,  0.98078528,
                            -0.70710677, -0.19509033,  0.92387954, -0.83146960},
                          {  1.00000000, -0.83146961,  0.38268343,  0.19509033,
                            -0.70710679,  0.98078528, -0.92387953,  0.55557022},
                          {  1.00000000, -0.98078528,  0.92387953, -0.83146961,
                             0.70710677, -0.55557022,  0.38268341, -0.19509030}
                        };

 fix12.10 costab2[:,:,:,:] =
      for elem1 in costab cross elem2 in costab
      return (array(elem1*elem2));

 // ***** Computes the Cosine transform for 8x8 windows *****
 $OutType result[:,:] =
      for window win[8,8] in src step (8,8)
           { $OutType wintransf[:,:] =
                  for VsipIndexesType u, VsipIndexesType v  in [8,8]
                       { float sumval =
                              for elem in win at (uint4 i, uint4 j)
                              return (sum(($InterType)(($InterType)elem *
                                           costab2[u,i,v,j] * alpha2[i,j])));

                         $OutType transf = AuxRoundInt(sumval);
                       }
                  return (array(transf));
           }
      return(tile(wintransf));

} return(result);

$end_generic
 
/****************************************************************************
 *
 *  Function    : Vsip8x8InvDiscreteCosineTransformFloat
 *
 *  Description : Compute the 8 x 8 Inverse Discrete Cosine Transform of an
 *                Image. The Image must have size multiple of 8x8
 *
 *  Parameters  : src - input image object
 *
 *  Note        : This routine differ from the bit/byte/int because it does not
 *                round the result.
 *                The IntType is computed using the IOtype and assuming
 *                the actual alpha values (N = 8)
 *
 ****************************************************************************/
$generic

   $Vsip8x8InvDiscreteCosineTransformFloat = 'Vsip8x8InvDiscreteCosineTransformFloat',
        $IOType = 'float', $InterType = 'float';

$in

$IOType [:,:] $Vsip8x8InvDiscreteCosineTransformFloat($IOType src[:,:])
{

  // ***** recover the src size *****
  VsipIndexesType rs, VsipIndexesType cs = extents(src);

 // ***** Test size of the input image (src) (see header) *****
 assert(((rs > 0) && (cs > 0)),
         "ERROR: Image can not have zero dimention. (",rs,"x",cs,")\n");
 assert((rs % 8 == 0) && (cs % 8 == 0),
        "ERROR: Image must be multiple of 8x8.(",rs,"x",cs,")\n");

 // ***** set alpha. alpha(0) = sqrt(1/N).          *****
 // ***** alpha(u) = sqrt(2/N). In this case N = 8  *****
 fix11.10 alpha[8]= {0.35355339, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5};

 fix11.10 alpha2[:,:] =
      for elem1 in alpha cross elem2 in alpha
      return (array(elem1*elem2));

 // ***** Table of precomputed values used inside the convolution      *****
 // ***** cos((2 * i + 1) * u * PI / (2 * N)) for i, u = 0, ..., 7     *****
 // ***** Computed using PI = 3.14159265 and N = 8                     *****
 fix12.10 costab[8,8] = { {  1.00000000,  0.98078528,  0.92387953,  0.83146961,
                             0.70710678,  0.55557023,  0.38268343,  0.19509032},
                          {  1.00000000,  0.83146961,  0.38268343, -0.19509032,
                            -0.70710678, -0.98078528, -0.92387953, -0.55557024},
                          {  1.00000000,  0.55557023, -0.38268343, -0.98078528,
                            -0.70710678,  0.19509032,  0.92387953,  0.83146962},
                          {  1.00000000,  0.19509032, -0.92387953, -0.55557024,
                             0.70710678,  0.83146962, -0.38268342, -0.98078528},
                          {  1.00000000, -0.19509032, -0.92387953,  0.55557023,
                             0.70710679, -0.83146961, -0.38268344,  0.98078528},
                          {  1.00000000, -0.55557023, -0.38268344,  0.98078528,
                            -0.70710677, -0.19509033,  0.92387954, -0.83146960},
                          {  1.00000000, -0.83146961,  0.38268343,  0.19509033,
                            -0.70710679,  0.98078528, -0.92387953,  0.55557022},
                          {  1.00000000, -0.98078528,  0.92387953, -0.83146961,
                             0.70710677, -0.55557022,  0.38268341, -0.19509030}
                        };

 fix12.10 costab2[:,:,:,:] =
      for elem1 in costab cross elem2 in costab
      return (array(elem1*elem2));

 // ***** Computes the Cosine transform for 8x8 windows *****
 $IOType result[:,:] =
      for window win[8,8] in src step (8,8)
           { $IOType wintransf[:,:] =
                  for VsipIndexesType u, VsipIndexesType v  in [8,8]
                       { float sumval =
                              for elem in win at (uint4 i, uint4 j)
                              return (sum(($InterType)(($InterType)elem *
                                           costab2[u,i,v,j] * alpha2[i,j])));

                         $IOType transf = sumval;
                       }
                  return (array(transf));
           }
      return(tile(wintransf));

} return(result);

$end_generic