//------------------------------------------------------------------------------
// Main data class		    																	--
// (C) Piero Giubilato 2008-2010, Berkeley Lab										   --
//------------------------------------------------------------------------------

//______________________________________________________________________________
// {Trace}
// [File name] "data_Layer.cpp"
// [Author] "Piero Giubilato"
// [Version] "1.1"
// [Modified by] "Piero Giubilato"
// [Last revision] "29 Jul 2009"
// [Language] "C++"
// [Compiler] "Visual C++ 9.x"
// [Member of] "Cool SEAL"
// [Project] "SEAL"
// [Description] "Main data class"
// [Key documentation]
// "Root user's guide"
// "Visual C++ Reference Help"
// {Trace}
//______________________________________________________________________________

// Overloading check
#ifndef data_Layer_C
#define data_Layer_C

// Standard components
#define _CRT_RAND_S
#include <stdlib.h>
#include <stdio.h>
#include <limits.h>
#include <time.h>
#include <math.h>

// Application components
#include "data_Layer.h"

// *****************************************************************************
// ****									  SPECIAL MEMBERS									 ****
// *****************************************************************************

//______________________________________________________________________________
template<class T>
data_Layer<T>::data_Layer(): data_Object(),
					dp_Array(0), dp_array_col_Count(0), dp_array_row_Count(0)
{
	// Provides the base class for raw data storage.
	// A layer can handle a 2D, 1D or 0D array of data of any arbitrary type.
	// Data addressing is implemented both via col, row indexes and via global
	// Idx indexing.
	// A layer object is a complete data type, i.e. it provides complete copy
	// constructor and assign operator; both defined for same type or type
	// casting operations.
   //
   // WARNING: this class provides also low-level data access methods (bare 
   // pointers) to save memory space and access time: this implies that memory
   // corruption errors may arise if the access is not correctly implemented!

	// Debug
	dbg_Print("data_Layer::data_Layer:()", DBG_LVL_ZERO);

   // Data preset with no parameters (1 data point)
	data_Resize(1, 1);
}

//______________________________________________________________________________
template<class T>
data_Layer<T>::data_Layer(UInt_t col_Count, UInt_t row_Count): data_Object(), 
					dp_Array(0), dp_array_col_Count(0), dp_array_row_Count(0)
{
	// Constructor with assigned number of cols and rows
	dbg_Print("data_Layer::data_Layer:(col, row)", DBG_LVL_ZERO);

   // Sizes the data array (set col and rows)
   data_Resize(col_Count, row_Count);
}

//______________________________________________________________________________
template<class T>
data_Layer<T>::data_Layer(const data_Layer<T>& layer_Ref): data_Object(), 
					dp_Array(0), dp_array_col_Count(0), dp_array_row_Count(0)
{
	// Default copy construnctor (to avoid compiler override)
	dbg_Print("data_Layer::data_Layer:(data_layer<T>&)", DBG_LVL_ZERO);

   // Sizes the data array (set col and rows)
   data_Copy(layer_Ref);
}

//______________________________________________________________________________
template<class T> template<class U>
data_Layer<T>::data_Layer(const data_Layer<U>& ref): data_Object(), 
					dp_Array(0), dp_array_col_Count(0), dp_array_row_Count(0)
{
	// Casting copy constructor
	dbg_Print("data_Layer::data_Layer:(U&)", DBG_LVL_ZERO);

   // Sizes the data array (set col and rows)
   data_Copy(ref);
}

//______________________________________________________________________________
template<class T>
data_Layer<T>::~data_Layer ()
{
   // Standard destructor
	dbg_Print("data_Layer::~data_Layer:()", DBG_LVL_ZERO);
	
	// Decreases the global memory allocation counter
	data_Object::stc_data_size_Del(dp_Size());
	
	// Frees all allocated memory
   delete []dp_Array;
}


// *****************************************************************************
// ****							    PRIVATE OBJECT ACCESS		 						 ****
// *****************************************************************************

//______________________________________________________________________________
template<class T>
void data_Layer<T>::data_Resize(UInt_t col_Count, UInt_t row_Count) 
{
   // Sets the size of the data container (dp_Array) accordingly to the given size.
	dbg_Print("data_Layer::data_Resize:(UInt, UInt)", DBG_LVL_STEP);	
	
   // Decreases the global memory allocation counter
	data_Object::stc_data_size_Del(sizeof(T) * dp_array_col_Count * dp_array_row_Count);

	// Check
   if (col_Count == 0) col_Count = 1;
   if (row_Count == 0) row_Count = 1;
	
	// Clears present mem locations
	if (dp_Array) delete []dp_Array;

	// Allocates the new array
	if (!(dp_Array = new T[col_Count * row_Count])) throw DBG_WHERE;
   
	// Increases the global memory allocation counter
	data_Object::stc_data_size_Add(sizeof(T) * col_Count * row_Count);
	
	// Stores the new dimensions
   dp_array_col_Count = col_Count;
   dp_array_row_Count = row_Count;
}

//______________________________________________________________________________
template<class T> template<class U>
void data_Layer<T>::data_Copy(const data_Layer<U>& Source)
{
	// Copies the data from 'Source', casting them. Does not copy anything
	// if the source is the same as the target.
	dbg_Print("data_Layer::data_Copy:(data_Layer<U>&)", DBG_LVL_STEP);

	// Self copy check
	if (this == (data_Layer<T>*)&Source) return;
	
	// Parent copy
	data_Object::operator=(Source);
	
	// Resizes the dp_Array only if necessary
	if((dp_Array == 0) || (Source.dp_col_Count() != dp_col_Count()) || (Source.dp_row_Count() != dp_row_Count())) 
		data_Resize(Source.dp_col_Count(), Source.dp_row_Count());
	
	// Deep and fast copy
	T* tF = this->dp_Front(); 
	U* sF = Source.dp_Front();
	UInt_t iCnt = Source.dp_Count(); 
	for (UInt_t i = 0; i < iCnt; i++) tF[i] = (T)sF[i];
}

//______________________________________________________________________________
template<class T> template<class U> 
void data_Layer<T>::data_Math(const U& Source, UInt_t Operation)
{
	// Performs the operation specified by 'Operation' by applying Scalar to
	// '*this', accordingly to Operation and applying cast.
	dbg_Print("data_Layer::data_Math:(U&, UInt)", DBG_LVL_STEP);

	// Checks for division by zero, in case returns a 0 filled layer
	if (Operation == km_Div && Source == 0) return;
	
	// Creates a target object for new assignement
	T* tF = this->dp_Front(); 
	UInt_t iCnt = dp_Count(); 
		
	// Performs the scalar algebra
	switch(Operation) {
		case km_Set: for (UInt_t i = 0; i < iCnt; i++) tF[i] = (T)Source; break;
		case km_Add: for (UInt_t i = 0; i < iCnt; i++) tF[i] += (T)Source; break;
		case km_Sub: for (UInt_t i = 0; i < iCnt; i++) tF[i] -= (T)Source; break;
		case km_Mul: for (UInt_t i = 0; i < iCnt; i++) tF[i] *= (T)Source; break;
		case km_Div: for (UInt_t i = 0; i < iCnt; i++) tF[i] /= (T)Source; break;
		case km_Sqr: for (UInt_t i = 0; i < iCnt; i++) tF[i] = (T)sqrt((double)tF[i]); break;
		case km_Pow: for (UInt_t i = 0; i < iCnt; i++) tF[i] = (T)pow((double)tF[i], (int)Source); break;
		case km_Mir: for (UInt_t i = 0; i < iCnt; i++) tF[i] = (T)(Source - tF[i]); break;
	}
}

//______________________________________________________________________________
template<class T> template<class U>
void data_Layer<T>::data_Math(const data_Layer<U>& Source, UInt_t Operation)
{
	// Performs the operation specified by 'Operation' value between 'Source' and
	// '*this', applying cast. If the dimensions do not match, only the overlapping
	// part of  '*this' will be affected by the manipulation, while the remaining
	// will retain the same value it had before.
	dbg_Print("data_Layer::data_Math:(data_Layer<U>&, UInt)", DBG_LVL_STEP);

	// Determines the overlap Region boundaries
	UInt_t cCnt; 
	UInt_t rCnt; 
	UInt_t tcCnt = dp_array_col_Count;
	UInt_t trCnt = dp_array_row_Count;
	UInt_t scCnt = Source.dp_col_Count();
	UInt_t srCnt = Source.dp_row_Count();
	(tcCnt <= scCnt) ? cCnt = tcCnt : cCnt = scCnt;
	(trCnt <= srCnt) ? rCnt = trCnt : rCnt = srCnt;
	
	// Gets the arrays fronts
	T* tF = this->dp_Front(); 
	U* sF = Source.dp_Front();
	UInt_t iCnt = dp_Count(); 

	// If dimensions do not match, goes for an overlap area algebra
	if((scCnt != tcCnt) || (srCnt != trCnt)) {
		switch(Operation) {
			case km_Add: for (UInt_t col = 0; col < cCnt; col++) for (UInt_t row = 0; row < rCnt; row++)
							tF[row * tcCnt + col] += (T)sF[row * scCnt + col]; break;
			case km_Sub: for (UInt_t col = 0; col < cCnt; col++) for (UInt_t row = 0; row < rCnt; row++)
							tF[row * tcCnt + col] -= (T)sF[row * scCnt + col];	break;
			case km_Mul: for (UInt_t col = 0; col < cCnt; col++) for (UInt_t row = 0; row < rCnt; row++)
							tF[row * tcCnt + col] *= (T)sF[row * scCnt + col];	break;
			case km_Div: for (UInt_t col = 0; col < cCnt; col++) for (UInt_t row = 0; row < rCnt; row++)
							tF[row * tcCnt + col] /= (T)sF[row * scCnt + col];	break;
		}
	// Otherwise, uses a faster linear matching
	} else {
		switch(Operation) {
			case km_Add: for (UInt_t i = 0; i < iCnt; i++) tF[i] += (T)sF[i]; break;		
			case km_Sub: for (UInt_t i = 0; i < iCnt; i++) tF[i] -= (T)sF[i]; break;		
			case km_Mul: for (UInt_t i = 0; i < iCnt; i++) tF[i] *= (T)sF[i]; break;		
			case km_Div: for (UInt_t i = 0; i < iCnt; i++) tF[i] /= (T)sF[i]; break;
		}
	}
}

//______________________________________________________________________________
template<class T>
void data_Layer<T>::dp_all_Stat(Double_t& avg, Double_t& sgm) const
{
	// Performs basic stat calculations
	dbg_Print("data_Layer::dp_all_Stat:()", DBG_LVL_FLOW);
	
	// Returns the minimum value present into the array
   UInt_t i_Count = dp_Count();
	//T min = dp_Array[0];
	//T max = dp_Array[0];
	Double_t sum = (Double_t)dp_Array[0];
	Double_t sum2 = pow((Double_t)dp_Array[0], 2);
	avg = (Double_t)dp_Array[0];
	for (UInt_t i = 1; i < i_Count; i++) {
		//if (dp_Array[i] < min) min = dp_Array[i];
		//if (dp_Array[i] > max) max = dp_Array[i];
		avg += (Double_t)dp_Array[i];
		sum += (Double_t)dp_Array[i];
		sum2 += pow((Double_t)dp_Array[i], 2);
	}

	// Calculates results
	avg /= i_Count;
	sgm = sqrt(pow(sum, 2) + (sum2 /= I_Count));
}

//______________________________________________________________________________
template<class T>
void data_Layer<T>::dp_all_Limits(T& min, T& max) const
{
	// Gets the max e min values inside the layer
	dbg_Print("data_Layer::dp_all_Limit:()", DBG_LVL_FLOW);
	
	// Returns the minimum value present into the array
   min = dp_Array[0];
	max = dp_Array[0];
	
	// Find the min, max
	UInt_t i_Count = dp_Count();
	for (UInt_t i = 1; i < i_Count; i++) {
		if (dp_Array[i] < min) min = dp_Array[i];
		if (dp_Array[i] > max) max = dp_Array[i];
	}
}

//______________________________________________________________________________
template<class T>
UInt_t data_Layer<T>::obj_Type() const
{
	// Returns the object type (the class type itself)
   return data_Object::kot_Layer; 
}

//______________________________________________________________________________
template<class T>
UInt_t data_Layer<T>::data_Type() const
{
	// Returns the data type (the object storage data type)
   dbg_Print("data_Layer::data_Type:()", DBG_LVL_STEP);

   // Try matching
   if(typeid(T) == typeid(UShort_t)) return data_Object::kdt_UShort;
   if(typeid(T) == typeid(Short_t)) return data_Object::kdt_Short;
   if(typeid(T) == typeid(UInt_t)) return data_Object::kdt_UInt;
   if(typeid(T) == typeid(Int_t)) return data_Object::kdt_Int;
   if(typeid(T) == typeid(Float_t)) return data_Object::kdt_Float;
   if(typeid(T) == typeid(Double_t)) return data_Object::kdt_Double;
	
	// Do we really need those?
	//if(typeid(T) == typeid(unsigned short)) return data_Object::kdt_Short;
	//if(typeid(T) == typeid(short)) return data_Object::kdt_Short;
	//if(typeid(T) == typeid(int)) return data_Object::kdt_Int;
	//if(typeid(T) == typeid(float)) return data_Object::kdt_Int;
	//if(typeid(T) == typeid(double)) return data_Object::kdt_Double;
	   
   // No matching data, returns null
   return data_Object::kdt_Null; 
}

// *****************************************************************************
// ****									  SIZE SET/READ		 							 ****
// *****************************************************************************

//______________________________________________________________________________
template<class T>
inline UInt_t data_Layer<T>::dp_col_Count() const
{
	// Returns the number of columns of px_Array.
	return dp_array_col_Count;
}

//______________________________________________________________________________
template<class T>
UInt_t data_Layer<T>::dp_col_Count(UInt_t col_Count)
{
	// Sets the number of columns of dp_Array. Attention: this will 
   // RESET the dp_Array! Returns the new number of columns, 0 on error.
	dbg_Print("data_Layer::dp_col_Count:(UInt)", DBG_LVL_STEP);

   // Check
   if (col_Count == 0) return dp_array_col_Count;

   // Resizes
   data_Resize(col_Count, dp_array_row_Count);
   
   // Return col number
   return dp_array_col_Count;
}

//______________________________________________________________________________
template<class T>
inline UInt_t data_Layer<T>::dp_row_Count() const
{
	// Returns the number of rows of px_Array.
	return dp_array_row_Count;
}

//______________________________________________________________________________
template<class T>
UInt_t data_Layer<T>::dp_row_Count(UInt_t row_Count)
{
	// Sets the number of rows and rows of dp_Array. Attention: this will 
   // RESET the dp_Array! Returns the new number of rows, 0 on error.
	dbg_Print("data_Layer::dp_row_Count:(UInt)", DBG_LVL_STEP);

   // Check
   if (row_Count == 0) return dp_array_row_Count;

   // Resizes
   data_Resize(dp_array_col_Count, row_Count);
   
   // Return row number
   return dp_array_row_Count;
}

//______________________________________________________________________________
template<class T>
inline UInt_t data_Layer<T>::dp_Count() const
{
	// Returns the current number of pixels stored.
   return  dp_array_col_Count * dp_array_row_Count;
}

//______________________________________________________________________________
template<class T>
UInt_t data_Layer<T>::dp_Count(UInt_t Count)
{
	// Sets the number of data point assuming a 1-dim array (1 row). This will 
   // RESET the dp_Array! Returns the new number of data points, 0 on error.
	dbg_Print("data_Layer::dp_Count:(UInt)", DBG_LVL_STEP);

   // Checks
   if (Count == 0) return dp_Count();
   
   // Resizes
   data_Resize(Count, 1);

	// Returns back the actual dp_Count
	return dp_Count();
}

//______________________________________________________________________________
template<class T>
UInt_t data_Layer<T>::dp_Count(UInt_t col_Count, UInt_t row_Count)
{
	// Sets the number of columns and rows of dp_Array. Attention: this will 
   // RESET the dp_Array! Returns the new number of data points, 0 on error.
	dbg_Print("data_Layer::dp_Count:(UInt, UInt)", DBG_LVL_STEP);

   // Checks
   if (col_Count == 0) return dp_Count();
   if (row_Count == 0) return dp_Count();

   // Resizes dp_Array
   data_Resize(col_Count, row_Count);
   
   // All OK
   return dp_Count();
}


// *****************************************************************************
// ****									    DATA ACCESS		 							 ****
// *****************************************************************************

//______________________________________________________________________________
template<class T>
inline UInt_t data_Layer<T>::dp_Idx(UInt_t Col, UInt_t Row) const
{
	// Returns the indexed (Idx) address of a data point from its
	// column and row address
	return Row * dp_array_col_Count + Col;
}

//______________________________________________________________________________
template<class T>
inline UInt_t data_Layer<T>::dp_Col(UInt_t Idx) const
{
	// Returns the column address of a data point from its index
	return (Idx % dp_array_col_Count);
}

//______________________________________________________________________________
template<class T>
inline UInt_t data_Layer<T>::dp_Row(UInt_t Idx) const
{
	// Returns the row address of a data point from its index
	return ((Idx - (Idx % dp_array_col_Count)) / dp_array_col_Count);
}

//______________________________________________________________________________
template<class T>
inline T data_Layer<T>::dp_Value(UInt_t Idx) const
{
	// Retrives the value of a data point at a given position, safe access
   // method. By using the direct memory access through the dp_Front 
   // reference value will give a 3 times faster data access.

   // Checks for out of bounds
   if (Idx >= dp_Count()) Idx = 0;
   
   // Returns the value
	return dp_Array[Idx];
}

//______________________________________________________________________________
template<class T>
inline T& data_Layer<T>::dp_Value(UInt_t Idx)
{
	// Checks for out of bounds
   if (Idx >= dp_Count()) Idx = 0;
   
   // Returns the value
	return dp_Array[Idx];
}

//______________________________________________________________________________
template<class T>
inline T data_Layer<T>::dp_Value(UInt_t Col, UInt_t Row) const
{
	// Retrives the value of a data point at a given position, safe access
   // method. By using the direct memory access through the dp_Front 
   // reference value will give a 4.5 times faster data access, at the price
   // of no safety check and col/row addressing.

   // Check for out of bounds
   if (Col >= dp_array_col_Count) Col = 0;
   if (Row >= dp_array_row_Count) Row = 0;

   // Returns the value
	return dp_Array[Row * dp_array_col_Count + Col];
}

//______________________________________________________________________________
template<class T>
inline T& data_Layer<T>::dp_Value(UInt_t Col, UInt_t Row)
{
	// Check for out of bounds
   if (Col >= dp_array_col_Count) Col = 0;
   if (Row >= dp_array_row_Count) Row = 0;

   // Returns the value
	return dp_Array[Row * dp_array_col_Count + Col];
}

//______________________________________________________________________________
template<class T> 
data_Layer<T>& data_Layer<T>::dp_all_Range(T value_Start, T value_Step)
{
	// Sets all the dp starting from 'value_Start' and increasing by 
   // 'value_Step' at each point.
   UInt_t i_Count = dp_Count();
   for (UInt_t i = 0; i < i_Count; i++) {
      dp_Array[i] = value_Start;
      value_Start += value_Step;
   }

	// Returns itself
	return *this;
}

//______________________________________________________________________________
template<class T> 
data_Layer<T>& data_Layer<T>::dp_all_Lattice(Int_t radius, Int_t space)
{
	// Fills the layer with a pseudo-lattice blob-like pattern
   UInt_t i_Count = dp_Count();
	UInt_t col_Count = dp_col_Count();
	UInt_t row_Count = dp_row_Count();
	Int_t pX = radius, pY = radius;
	Double_t r = 0;

	// Cleans
	for (UInt_t i = 0; i < i_Count; i++) dp_Array[i] = 0;
	
	// Fill
	while (pY < row_Count - radius) {
		while (pX < col_Count - radius) {
			for (Int_t y = pY - radius; y < pY + radius; y++) {
				for (Int_t x = pX - radius; x < pX + radius; x++) {
					r = sqrt((Double_t)(pX - x)*(pX - x) + (pY - y)*(pY - y));
					if (r < radius) dp_Value(x, y) = (T)(255 * (1 - r/radius));
				}
			}
			pX += space;
		}
		pX = radius;
		pY += space;
	}

	// Returns itself
	return *this;
}

//______________________________________________________________________________
template<class T> 
data_Layer<T>& data_Layer<T>::dp_all_Wave(Double_t lenght, Double_t tilt)
{
	// Fills the layer with a sine-wave of wavelenght lenght and tilt oriented
   UInt_t i_Count = dp_Count();
	UInt_t col_Count = dp_col_Count();
	UInt_t row_Count = dp_row_Count();
	
	// Cleans
	for (UInt_t i = 0; i < i_Count; i++) dp_Array[i] = 0;
	
	// Fill
	for (Int_t y = 0; y < row_Count; y++) {
		for (Int_t x = 0; x < col_Count; x++) {
			dp_Array[x + y * col_Count] = (T)(128 + 128 * sin(x * lenght));
		}
	}

	// Returns itself
	return *this;
}

//______________________________________________________________________________
template<class T> 
data_Layer<T>& data_Layer<T>::dp_all_Populate(UInt_t cl_Count, UInt_t cl_Height, UInt_t cl_Size)
{
	// Populates the layer with clusters
   // cl_Count: number of cluster
	// cl_Height: cluster average height
	// cl_Size: cluster average diameter
	
	UInt_t i_Count = dp_Count();
	UInt_t col_Count = dp_col_Count();
	UInt_t row_Count = dp_row_Count();
	Float_t x, y, radius, height, var;
	Float_t dst;
   	
	// Seeds random generator
	srand((unsigned)time(NULL));

   // Cleans
	for (UInt_t i = 0; i < i_Count; i++) dp_Array[i] = 0;
	
	// Populates
	for (UInt_t cl_Idx = 0; cl_Idx < cl_Count; cl_Idx++) {
		
		// Get cluster properties
		x = (Float_t)((double)rand() / (RAND_MAX + 1) * col_Count);
		y = (Float_t)((double)rand() / (RAND_MAX + 1) * row_Count);
		radius = (Float_t)((double)rand() / (RAND_MAX + 1) * cl_Size);
		height = (Float_t)((double)rand() / (RAND_MAX + 1) * cl_Height);	
		
		// Populates pixels
		for (Float_t b = - radius; b < radius; b++) {
			for (Float_t a = -radius; a < radius; a++) {
				if ((x + a) >= 0 && (x + a) < col_Count && (y + b) >= 0 && (y + b) < row_Count) {
					dst = (a*a + b*b);
					if (dst <= radius*radius) {
						var = (Float_t)((double)rand() / (RAND_MAX + 1) * cl_Height / 4);	
						dp_Value((UInt_t)(x + a), (UInt_t)(y + b)) = (T)(height - sqrt(dst)/radius*height + var);
					}
				}
			}
		}
	}
			
	// Returns itself
	return *this;
}
//______________________________________________________________________________
template<class T> 
data_Stack<T>* data_Layer<T>::dp_all_Clusterize(T min_Thr, T seed_Thr)
{
	// Creates a stack of clusters generated by the layer data. min_Thr is the 
	// minimum value to consider a dp, seed_Thr is the minimum value the highest
	// dp into a cluster must reach to consider the cluster
   
	// Dimensions
	UInt_t i_Count = dp_Count();
	UInt_t col_Count = dp_col_Count();
	UInt_t row_Count = dp_row_Count();
	UInt_t row_Idx = 0;
	UInt_t cl_Idx = 0;
	bool cl_On = false;
	bool cl_Close = false;
	UInt_t cl_Sx = 0;
	UInt_t cl_Ex = 0;
	
	// Pivot scanline array
	UInt_t* l_Array = new UInt_t[col_Count];
	for (UInt_t i = 0; i < col_Count; i++) l_Array[i] = 0;
	
	// The stack to be returned
	data_Stack<T>* stack = new data_Stack<T>;
	data_Point<T>* point = new data_Point<T>(0,0,0);

	// For each row
	for (UInt_t y = 0; y < row_Count; y++) {
		
		// Pivot indexes
		row_Idx = y * col_Count;
				
		// Scanline cluster builder
		for (UInt_t x = 0; x < col_Count; x++) {
			
			// Starts a cluster line
			if (dp_Array[row_Idx + x] > min_Thr) {
				
				// Starts the cluster line in case
				if (!cl_On) {
					cl_On = true;
					cl_Idx = 0;
					cl_Sx = x;
					
					// Check for the bl corner if possible
					if (x > 0) {if (l_Array[x - 1]) cl_Idx = l_Array[x - 1];}
					
					// Fast close in case
					if (x == col_Count - 1) cl_Close = true;

				// Invoke cluster close on last column
				} else if (x == col_Count - 1) cl_Close = true;
											
				// Check for any already existent neighborhood cluster
				if (l_Array[x]) {
					if (cl_Idx && (cl_Idx != l_Array[x])) stack->Link(cl_Idx - 1, l_Array[x] - 1);
					cl_Idx = l_Array[x]; // Join the cluster
				}
				
				// Sets the last good pixel position
				cl_Ex = x;
			
			// Close cluster line otherwise
			} else if (cl_On) cl_Close = true;
										
			// Close cluster line
			if (cl_Close) {
				
				// Close the scanline
				cl_Close = false;
				cl_On = false;
					
				// Grabs the last possible neighborod in case (br corner)
				if (cl_Ex < (col_Count - 1)) {
					if(l_Array[cl_Ex + 1]) {
						if (cl_Idx && (cl_Idx != l_Array[cl_Ex + 1])) stack->Link(cl_Idx - 1, l_Array[cl_Ex + 1] - 1);
						cl_Idx = l_Array[cl_Ex + 1];
					}
				}
					
				// Creates a new cluster if no one neighborhood found
				if (!cl_Idx) {
					stack->cl_Array.push_back(new data_Cluster<T>());
					cl_Idx = stack->cl_Array.size();
				}

				// Fills the cluster with the cluster line points
				for (UInt_t p = cl_Sx; p <= cl_Ex; p++) {
					stack->cl_Array[cl_Idx - 1]->dp_Array.push_back(new data_Point<T>(p, y, dp_Array[row_Idx + p]));
					l_Array[p] = cl_Idx;
				}

				// Clears the scanline only if the end pixel
				// is different from the curren tone
				if (cl_Ex != x) l_Array[x] = 0;

			} else l_Array[x] = 0; // Clears the scanline
		}
	}

	// Removes the pivot linescan array
	delete[] l_Array;
	
	// Connects linked clusters
	stack->Clean();

	// Returns the fulfilled stack
	return stack;
}

//______________________________________________________________________________
template<class T> inline T* data_Layer<T>::dp_Front() const
{
	// Returns the address of dp_Array.
	return dp_Array;
}

//______________________________________________________________________________
template<class T> size_t data_Layer<T>::dp_Size() const
{
	// Returns the current size of dp_Array (bytes).
	return  dp_Count() * sizeof(T);
}

// *****************************************************************************
// ****								    ROOT OVERLOADS		 							 ****
// *****************************************************************************

//______________________________________________________________________________
template<class T> void data_Layer<T>::Dump(UInt_t level) const
{
   // Dumps all the object data into std::out
   TString space(' ', level);
   
   // Shows the header 
	std::cout << space.Data() << "*data_Layer: " << (int*)this << "\n";
   std::cout << space.Data() << "class_Ver: " << class_Ver << "\n";
   
	// Show the parent object
   data_Object::Dump(level + 2);
	
	// Shows and the data members
   std::cout << space.Data() << "col_Count: " << dp_array_col_Count << "\n";
   std::cout << space.Data() << "row_Count: " << dp_array_row_Count << "\n";

   // Show the data matrix member
   for (UInt_t row = 0; row < dp_array_row_Count; row++) {
      std::cout << space.Data();
      for (UInt_t col = 0; col < dp_array_col_Count; col++) {
         std::cout << dp_Value(col, row) << " ";
		}
      std::cout << "\n";
   }
}

//______________________________________________________________________________
template<class T> void data_Layer<T>::Streamer(TBuffer& b)
{
   // Fills the standar root streamer with the current data points

   // Read
   if (b.IsReading()) {
      
      // Debug
      dbg_Print ("data_Layer::Streamer:(read TBuffer&)", DBG_LVL_STEP)

      // Reads class version
      UShort_t class_ver_Read = 0;
      b >> class_ver_Read;
         
      // Check version
      if (class_ver_Read != class_Ver) {
         
         // Debug
         dbg_Print ("data_Layer::Streamer:(read): version mismatch", DBG_LVL_WARN)
         dbg_Value(class_ver_Read, DBG_LVL_WARN);
         dbg_Value(class_Ver, DBG_LVL_WARN);
      }

      // Reads parent class
      data_Object::Streamer(b); 

      // Reads data points
      b >> dp_array_col_Count;
      b >> dp_array_row_Count;

      // Allocates the needed storage space
      data_Resize(dp_array_row_Count, dp_array_col_Count);
      
      // Read the data into the matrix
      b.ReadFastArray(dp_Array, dp_array_row_Count * dp_array_col_Count);

   // Writes
   } else {
      
      // Debug
      dbg_Print ("data_Layer::Streamer:(write TBuffer&)",DBG_LVL_STEP)

      // Writes class version
      b << class_Ver;

      // Writes parent class
      data_Object::Streamer(b);

      // Writes the data points
      b << dp_array_col_Count;
      b << dp_array_row_Count;
      b.WriteFastArray(dp_Array, dp_array_row_Count * dp_array_col_Count);
   }
}

// Overloading check
#endif