/*********************************************************************
 * This file is part of the PRAPI library.
 *
 * Copyright (C) 2001 Topi Mäenpää
 * All rights reserved.
 *
 * This program is free software. You can redistribute and/or modify
 * it under the terms of the free software licence found in the
 * accompanying file "COPYING". The licence terms must always be
 * redistributed with this source file. The above copyright notice
 * must be reproduced in all modified and unmodified copies of this
 * source file.
 *
 * $Revision$
 *********************************************************************/

#include "PCExtractor.h"

#include "../dsp/Wavelet.h"

using namespace util;
using namespace prapi::dsp;
	
namespace prapi { namespace extras {

	List<double> PCExtractor::getFeatureVector(const Matrix<double>& mat)
		throw (FeatureExtractionException&)
	{
		int r = mat.getRows(), c = mat.getColumns();
		if (!_bSplit || (r == 16 && c == 16))
			{
				List<Matrix<double> > decomp1(4);
				List<Matrix<double> > decomp2(4);
				List<Matrix<double> > decomp3(4);
				
				Wavelet::dwt<double>(mat,"db2",decomp1);
				Wavelet::dwt(decomp1[0],"db2",decomp2);
				Wavelet::dwt(decomp2[0],"db2",decomp3);
				
				List<double> result(4);
				result += l1norm(decomp3[0]);
				result += l1norm(decomp3[1]) + l1norm(decomp3[2]);
				result += l1norm(decomp2[1]) + l1norm(decomp2[2]);
				result += l1norm(decomp1[1]) + l1norm(decomp1[2]);

				return result;
			}
		else
			{
				int nx = c>>4, ny = r>>4, ox = (c-(nx<<4))>>1, oy = (r-(ny<<4))>>1;
				List<Matrix<double> > pieces(MatrixUtils::split(mat,16,16,oy,ox));
				List<double> result(4);
				result.setLength(4);
				result = 0.0;
				for (int i=0;i<pieces.getLength();i++)
					result.add(getFeatureVector(pieces[i]));
				result /= pieces.getLength();

				return result;
			}

	}

	double PCExtractor::l1norm(const Matrix<double>& mat)
	{
		double sum = 0;
		const double* data = mat.getData();
		int size = mat.getRows() * mat.getColumns();
		for (int i=size;i--;data++)
			sum += absolute(*data);
		return sum/size;
	}
	
}}