/*
Implementation of CCSDS 122.0-B-1 Recommended Standard
Please note:
(1)	Before you download and use the program, you must read and agree the license agreement carefully. 
(2)	We supply the source code and program WITHOUT ANY WARRANTIES. The users will be responsible 
        for any loses or damages caused by the use of the source code and the program. 

Author: 
Hongqiang Wang
Department of Electrical Engineering
University of Nebraska-Lincoln
Email: hqwang@bigred.unl.edu, hqwang@eecomm.unl.edu

Your comment and suggestions are welcome. Please report bugs to me via email and I would greatly appreciate it. 
Jan. 21, 2007
*/ 


#include <malloc.h>
#include <math.h>
#include "global.h"

extern void BlockScanEncode(StructCodingPara *PtrCoding, 	 BitPlaneBits *BlockInfo);

extern void StagesEnCoding(StructCodingPara *PtrCoding,   BitPlaneBits *BlockInfo);

extern void StagesDeCoding(StructCodingPara *PtrCoding,   BitPlaneBits *BlockInfo);

void ACDepthEncoder(StructCodingPara *PtrCoding,	BitPlaneBits *BlockInfo);

short ACDepthDecoder(StructCodingPara *PtrCoding,	 BitPlaneBits *BlockInfo);

void ACGaggleEncoding(StructCodingPara *PtrCoding, 
					  BitPlaneBits *BlockInfo, 
					  int StartIndex,
					  int gaggles, 
					  int Max_k, 
					  int ID_Length)
{
	int i = 0; 
	int k = 0; 
	int min_k = 0;
	int temp_DC = 0;
		
	DWORD32 min_bits = 0xFFFF; 
	DWORD32 total_bits = 0; 

	if (gaggles == 0)
		return;
	
	if(PtrCoding->PtrHeader->Header.Part3.OptACSelect == TRUE)
	{
		//calculate  the uncoded:
		WORD16 uncoded_bits = PtrCoding->N * gaggles;
		for ( k = 0; k <= Max_k; k ++)
		{
			total_bits = k * gaggles;

			for (i = StartIndex; i < StartIndex + gaggles; i ++)
				total_bits += ((BlockInfo[i].MappedAC >> k ) + 1);			
			
			if( total_bits < min_bits)
			{
				min_bits = total_bits;
				min_k = k;
			}
		}

		if (min_bits < uncoded_bits)
		{
			// Get the one that gives the minimum number of bits. 
			BitsOutput(PtrCoding, min_k, ID_Length);
			if(StartIndex == 1)
				BitsOutput(PtrCoding, BlockInfo[0].MappedAC, PtrCoding->N);
			for (i = StartIndex; i < StartIndex + gaggles; i ++)
			{
				temp_DC = BlockInfo[i].MappedAC;
				temp_DC >>= min_k;	
				temp_DC++;
				// determine the length of the symbol
				BitsOutput(PtrCoding, 1, temp_DC);	
			}
			
			for (i = StartIndex; i < StartIndex + gaggles; i ++)
				BitsOutput(PtrCoding, BlockInfo[i].MappedAC, min_k);
		}
		else
		{		
			UCHAR8 ones = ~0;			
			BitsOutput(PtrCoding, ones, ID_Length);			
			if(StartIndex == 1)
				BitsOutput(PtrCoding, BlockInfo[0].MappedAC, PtrCoding->N);
			for (i = StartIndex; i < StartIndex + gaggles; i ++)
				BitsOutput(PtrCoding, BlockInfo[i].MappedAC, PtrCoding->N);
			return;
		}
	}
	else
	{
		// Using so called Heuristic procedure. 
		int delta = 0;
		for (i = StartIndex; i < StartIndex + gaggles; i ++)
			delta += BlockInfo[i].MappedAC;
		if (64 *  delta >= 23 * gaggles * PtrCoding->N)
		{
			// uncoded/
			UCHAR8 ones = ~0;			
			BitsOutput(PtrCoding, ones, ID_Length);
			for (i = StartIndex; i < StartIndex + gaggles; i ++)
				BitsOutput(PtrCoding, BlockInfo[i].MappedAC, PtrCoding->N);
			return;
		}
		else if (207 * gaggles > 128 * delta)
		{
			k = 0;
		}
		else if (gaggles * (1<<(PtrCoding->N - 5)) <= 128 * delta + 49 * gaggles)
		{
			k = PtrCoding->N - 2;
		}
		else
		{
			long sum = 128 * delta + 49 * gaggles; 
			k = 0;
			while (gaggles * (1 << (k + 7)) <= sum)
				k++;
			k--;
		}	
		min_k = k; 
		BitsOutput(PtrCoding, min_k, ID_Length);

		if(StartIndex == 1)
			BitsOutput(PtrCoding, BlockInfo[0].MappedAC, PtrCoding->N);

		for (i = StartIndex; i < StartIndex + gaggles; i ++)
		{
			temp_DC = BlockInfo[i].MappedAC;
			temp_DC >>= min_k;	
			temp_DC++;
			// determine the length of the symbol
			BitsOutput(PtrCoding, 1, temp_DC);	
		}	
		for (i = StartIndex; i < StartIndex + gaggles; i ++)
			BitsOutput(PtrCoding, BlockInfo[i].MappedAC, min_k);
	}	
}




void DPCM_ACMapper(BitPlaneBits *BlockInfo,
				   int size,
				   short N)
{
	short *diff_AC;
	short theta = 0;  
	long i;
	int X_Min = 0;
	int X_Max = ((1 << N) - 1) ;

	diff_AC = (short *)calloc(size,sizeof(short)); 
	diff_AC[0] = BlockInfo[0].BitMaxAC;
	BlockInfo[0].MappedAC = BlockInfo[0].BitMaxAC;
	for(i = 1; i < size; i ++)
		diff_AC[i] = BlockInfo[i].BitMaxAC - BlockInfo[i-1].BitMaxAC;

	for ( i = 1; i < size; i ++)
	{
		theta = min(BlockInfo[i-1].BitMaxAC - X_Min, X_Max - BlockInfo[i-1].BitMaxAC);
		if (diff_AC[i] >= 0 && diff_AC[i] <= theta)
		{
			BlockInfo[i].MappedAC = 2 * diff_AC[i];
		}
		else if(diff_AC[i] < 0 && diff_AC[i] >= -theta)
		{
			BlockInfo[i].MappedAC = - 2 * diff_AC[i] - 1;
		}
		else
		{
			BlockInfo[i].MappedAC = theta + abs(diff_AC[i]);
		}
	}
	free(diff_AC);
	return;
}



void ACDepthEncoder(StructCodingPara *PtrCoding,
					BitPlaneBits *BlockInfo)
{
	SINT GaggleStartIndex; 
	UCHAR8 Max_k; 
	UCHAR8 ID_Length;
	UINT32 index = 0;
	UINT32 bit_num= 0;
	UINT32 counter = 0;
	UCHAR8 gaggles ; // initial value is 15.

	DOUBLE4 temp; 	

	temp = log10(1 + PtrCoding->PtrHeader->Header.Part1.BitDepthAC_5Bits) / log10(2);

	if (temp - (short)temp > 0 )			
	{
		PtrCoding->N = (UCHAR8)(temp + 1);
	}
	else
	{
		PtrCoding->N = (UCHAR8)temp;
	}

	DPCM_ACMapper(BlockInfo, PtrCoding->PtrHeader->Header.Part3.S_20Bits, PtrCoding->N);

	if ( PtrCoding->N == 2) 
	{
		Max_k = 0;
		ID_Length = 1;
	}
	else if (PtrCoding->N <=4)
	{
		Max_k = 2;
		ID_Length = 2;	
	}
	else if (PtrCoding->N <=5)
	{
		Max_k = 6;
		ID_Length = 3;	
	}
	else 
	{
		ErrorMsg(BPE_DATA_ERROR);	
	}

	// determine which one results in the minimum bits used. 
	// adaptive entropy coder. We will select the one results in
	// the minimum number of bits. 
	
	gaggles = GAGGLE_SIZE - 1;

	if (PtrCoding->PtrHeader->Header.Part3.S_20Bits < gaggles)
		gaggles = (UCHAR8) (PtrCoding->PtrHeader->Header.Part3.S_20Bits - 1);
	// output the first reference symbol:

	GaggleStartIndex = 1;

	ACGaggleEncoding(PtrCoding, BlockInfo, GaggleStartIndex,
		gaggles, Max_k, ID_Length);

	if(PtrCoding->SegmentFull == TRUE) 
		return;

	GaggleStartIndex += gaggles;
	gaggles = GAGGLE_SIZE;	

	while(PtrCoding->PtrHeader->Header.Part3.S_20Bits - GaggleStartIndex >= gaggles)
	{		 
		ACGaggleEncoding(PtrCoding, BlockInfo, GaggleStartIndex,	gaggles, Max_k, ID_Length);

		GaggleStartIndex += gaggles;

		if(PtrCoding->SegmentFull == TRUE) 
			return;
	}	
	gaggles = (UCHAR8)(PtrCoding->PtrHeader->Header.Part3.S_20Bits - GaggleStartIndex);		

	ACGaggleEncoding(PtrCoding, BlockInfo, GaggleStartIndex, gaggles, Max_k, ID_Length);

	return;
}



void ACBpeEncoding(StructCodingPara *PtrCoding,    BitPlaneBits *BlockInfo)
{
	WORD16 i = 0;
	UCHAR8 BitPlane = 0;

////////////////////////////////////////////////////////////////////////////////////////
	// The link nodes PlaneSym are used to coding of the bit planes. Contains a 
	// Link data structure. It encodes more than 16 blocks at this time. 
	// Need to modify it to encode up to 16 blocks only. May 20 2004
	if (PtrCoding->PtrHeader->Header.Part1.BitDepthAC_5Bits != 0) 
		// Bit plane scan is not necessary if BitDepthAC = 0. 
	{
		if (PtrCoding->PtrHeader->Header.Part1.BitDepthAC_5Bits == 1)
		{ //only the lowest bit plane has valud AC component manitude.
			for ( i = 0; i < PtrCoding->PtrHeader->Header.Part3.S_20Bits; i++)
				BitsOutput(PtrCoding, BlockInfo[i].BitMaxAC, 1);	
		}
		else
		{				
			// returned is the mapped ac value. 
			ACDepthEncoder(PtrCoding, BlockInfo);
		}
			
			// 3.2 Family tree scaning sequence in a segment. 

		for ( BitPlane = PtrCoding->PtrHeader->Header.Part1.BitDepthAC_5Bits ; 
		BitPlane > 0; BitPlane --)
		{
			PtrCoding->BitPlane = BitPlane;
			if((PtrCoding->PtrHeader->Header.Part2.BitPlaneStop_5Bits == BitPlane)
				&&(PtrCoding->PtrHeader->Header.Part1.Part2Flag == TRUE))
				return;

/////////////////////////////////////////////////////////////////////////
//////////////////////////  Stage 0   ///////////////////////////////////	
// encode the DC component (single bit only)
			if (PtrCoding->PtrHeader->Header.Part4.DWTType == INTEGER_WAVELET)
			{
				if ((BitPlane <= PtrCoding->QuantizationFactorQ) && 
					(PtrCoding->QuantizationFactorQ >
					PtrCoding->PtrHeader->Header.Part4.CustomWtLL3_2bits ) && 
					PtrCoding->PtrHeader->Header.Part4.CustomWtLL3_2bits < BitPlane)
				{						
					// output the DC bits that have been shifted out
					for(i = 0; i < PtrCoding->PtrHeader->Header.Part3.S_20Bits;i++)
						BitsOutput(PtrCoding, ((BlockInfo[i].DCRemainder >> (BitPlane - 1)) & 0x01), 1);
				}
			}
			else
			{
				if (BitPlane <= PtrCoding->QuantizationFactorQ) 
				{
					// output the DC bits that have been shifted out
					for(i = 0; i < PtrCoding->PtrHeader->Header.Part3.S_20Bits;i++)
						BitsOutput(PtrCoding, ((BlockInfo[i].DCRemainder >> (BitPlane - 1)) & 0x01), 1);
				}
			}
					
/////////////////////////////////////////////////////////////////////////					
		
			if(PtrCoding->SegmentFull == TRUE) 
				return;

			BlockScanEncode(PtrCoding, BlockInfo);	

			if(PtrCoding->SegmentFull == TRUE) 
				return;

			// ***************************  4. Stages coding *********************************//
			StagesEnCoding(PtrCoding, BlockInfo);		

			if(PtrCoding->SegmentFull == TRUE) 
				return;

		}		
	}
}



void CheckUsefill(StructCodingPara * PtrCoding)
{

	if(PtrCoding->PtrHeader->Header.Part2.SegByteLimit_27Bits != 0 
		&& PtrCoding->PtrHeader->Header.Part2.UseFill == TRUE)
	{
		if((PtrCoding->PtrHeader->Header.Part2.SegByteLimit_27Bits << 3) < 
			PtrCoding->Bits->SegBitCounter)
		{ 
			int RemainderBits = 0;