#ifndef IP_SYM_CONVOLVE
#define IP_SYM_CONVOLVE 1
////////////////////////////////////////////////////////////////////////////
//! docentry="Image.Image Processing.Spatial_Filters"
//! example=exIPConv.cc
//! file="amma/Image/IPStream/IPLinear/IPConvSym.hh"
//! lib=IPLinear
//! rcsid="$Id: IPConvSym.hh,v 1.1 2000/10/18 10:13:02 ees1wc Exp $"
//! author="Bill Christmas"
//! date="01/12/99"

#include "amma/DP/Process.hh"
#include "amma/NumImage.hh"
#include "amma/Image/IPSpatFilter.hh"
#include "amma/SArray1d.hh"

using namespace IPSpatFilterN;

// --------------------------------------------------------------------------
// **********  IPConvSymBodyC  ********************************************
// --------------------------------------------------------------------------

//! userlevel=Develop
//: Body class for IPConvSymC.

template <class PixelT, class CoeffT>
class IPConvSymBodyC
  : public DPProcessBodyC<NumImageC<PixelT>,NumImageC<PixelT> >,
    private IPSpatFilterC  {

public:
  IPConvSymBodyC (const SArray1dC<CoeffT> & Coeffs, DirnT Dirn,
		   ResizeT Resize, PadT Pad);
  //: Constructor
  // (See handle class IPConvSymC)

  virtual NumImageC<PixelT> Apply (const NumImageC<PixelT> &im);
  //:  Apply convolution to "im"
  // Uses mask specified in constructor.

protected:
  const SArray1dC<CoeffT> mask;  // filter coefficients
};



// --------------------------------------------------------------------------
// **********  IPConvSymC  ********************************************
// --------------------------------------------------------------------------

//! userlevel=Normal
//: Convolves an image with a 1-D odd-order symmetrical mask
//
// <p>Similar to <a href="IPConvolveC.html">IPConvolveC</a>, but with these
// differences:
//
// <ul>
// <li>The centre of the mask is used to align the input and output images
// <li>It is faster (hopefully)
// </ul>
//
// <p>Convolution can be vertical, horizontal or both. Because the mask is
// symmetrical, you only have to specify the right-hand &frac12; of it: the
// <b>central</b> coefficient of the mask is element 0 of the array.
//
// <p>The mask type is separately templated, so that you can use any
// combination of pixel and mask types, provided multiplication between them is
// supported.
//
// <p>See namespace <a href="IPSpatFilterN.html">IPSpatFilterN</a> for enumerated
// parameter values.

template <class PixelT, class CoeffT>
class IPConvSymC
  : public DPProcessC<NumImageC<PixelT>,NumImageC<PixelT> >
{
public:
  IPConvSymC();
  //: Default constructor.
  // Creates an invalid handle.
  
  IPConvSymC (const SArray1dC<CoeffT> & Coeffs, DirnT Direction,
	       ResizeT Resize=Erode, PadT Pad = PBRescale);
  //: Constructor with all params
  //
  // See IPSpatFilterN for explanation of parameter values
  //
  //!param: Coeffs - array of filter coefficients
  //!param: Otherwise - see IPSpatFilterC
};


/////////////////////////////////////////////////////////////////////////
// Implementation:

template <class PixelT, class CoeffT>
IPConvSymBodyC<PixelT,CoeffT>::IPConvSymBodyC (const SArray1dC<CoeffT> & Coeffs, DirnT Dirn, ResizeT Resize, PadT Pad)
  : IPSpatFilterC(((Dirn==Horizontal) ? 1 : 2*Coeffs.Size()-1),
		  ((Dirn==Vertical)   ? 1 : 2*Coeffs.Size()-1),
		  Dirn, Resize, Pad, ShiftUL),
    mask(Coeffs)
{}



template <class PixelT, class CoeffT>
NumImageC<PixelT> IPConvSymBodyC<PixelT,CoeffT>::Apply (const NumImageC<PixelT> &im)
{
  if (mask.Size() < 1) {
    errAMMA << "bad mask size: " << mask.Size();
    errAMMA.Function("NumImageC<PixelT> IPConvSymBodyC<PixelT,CoeffT>::Apply (const NumImageC<PixelT> &im)").Exit(-1);
  }
  if ((pad_method != PadBlack) && (pad_method != PBRescale)) {
    errAMMA << "Pad method " << pad_method << " not yet implemented";
    errAMMA.Function("NumImageC<PixelT> IPConvSymBodyC<PixelT,CoeffT>::Apply (const NumImageC<PixelT> &im)").Exit(-1);
  }
  NumImageC<PixelT> current;
  ImageRectangleC im_rect(im.Rectangle());
  PixelT Zero;  // create general-purpose zero value for mask coeff template
  StdMath::SetToZero(Zero);
  switch (resize) {
  case Erode:
    current = NumImageC<PixelT>(im);
    break;
  case Original:
    current = NumImageC<PixelT>(im_rect.Expand((rows-1)/2, (rows-1)/2,
					       (cols-1)/2, (cols-1)/2));
    current.Fill(Zero).Copy(im);
    break;
  case Expand:
    current = NumImageC<PixelT>(im_rect.Expand(rows-1, rows-1,
					       cols-1, cols-1));
    current.Fill(Zero).Copy(im);
    break;
  }
  if ((direction == Horizontal) || (direction == TwoD)) {
    NumImageC<PixelT> op(ImageRectangleC(current.TRow(), current.BRow(),
					 current.LCol()+mask.Size()-1,
					 current.RCol()-mask.Size()+1));
    if (op.IsValid()) {
      FOREACH_PIXEL (op, p) {
	op[p] = current[p] * mask[0];
	for (IndexC i=1; i<mask.Size(); ++i)
	  op[p] += (current[p.Row()][p.Col()-i]
		    + current[p.Row()][p.Col()+i]) * mask[i];
      }
      // if required, adjust weight of border pixels to preserve constant o/p
      // for constant i/p
      if (pad_method == PBRescale) {
	CoeffT tot(mask[0]); // sum of mask coeffs
	for (IndexC i=1; i<mask.Size(); ++i)  tot += 2*mask[i];
	CoeffT norm(tot); // sum of those coeffs still in o/p image
	for (IntT c = mask.Size()-2; im.LCol()+c >= op.LCol(); --c) {
	  norm-= mask[abs(c+1)]; // one more coeff falls outside o/p image
	  if (c < (IntT) im.Cnum()/2) {
	    FOREACH_ROW(op, r) op[r][im.LCol()+c] *= tot/norm;
	    FOREACH_ROW(op, r) op[r][im.RCol()-c] *= tot/norm;
	  }
	}
      }
      current = op;
    }
  }
  if ((direction == Vertical) || (direction == TwoD)) {
    NumImageC<PixelT> op(ImageRectangleC(current.TRow()+mask.Size()-1,
					current.BRow()-mask.Size()+1,
					current.LCol(), current.RCol()));
    if (op.IsValid()) {
      FOREACH_PIXEL (op, p) { 
	op[p] = current[p] * mask[0];
	for (IndexC i=1; i<mask.Size(); ++i)
	  op[p] += (current[p.Row()-i][p.Col()]
		    + current[p.Row()+i][p.Col()]) * mask[i];
      }
      // if required, adjust weight of border pixels to preserve constant o/p
      // for constant i/p
      if (pad_method == PBRescale) {
	CoeffT tot(mask[0]); // sum of mask coeffs
	for (IndexC i=1; i<mask.Size(); ++i)  tot += 2*mask[i];
	CoeffT norm(tot); // sum of those coeffs still in o/p image
	for (IntT r = mask.Size()-2; im.TRow()+r >= op.TRow(); --r) {
	  norm-= mask[abs(r+1)]; // one more coeff falls outside o/p image
	  if (r < (IntT) im.Rnum()/2) {
	    FOREACH_COL(op, c) op[im.TRow()+r][c] *= tot/norm;
	    FOREACH_COL(op, c) op[im.BRow()-r][c] *= tot/norm;
	  }
	}
      }
      current = op;
    }
  }
  return current;
}

template <class PixelT, class CoeffT>
IPConvSymC<PixelT,CoeffT>::IPConvSymC()
{}

template <class PixelT, class CoeffT>
IPConvSymC<PixelT,CoeffT>::
IPConvSymC (const SArray1dC<CoeffT> & Coeffs, DirnT Direction,
	     ResizeT Resize, PadT Pad)
  : DPProcessC<NumImageC<PixelT>,NumImageC<PixelT> >
(*new IPConvSymBodyC<PixelT,CoeffT>(Coeffs, Direction, Resize, Pad))
{}


#endif