spdec.h

GPS programme for changchun at 8X

//“星载大压缩比数据压缩技术”压缩编解码核心算法
// 预扫描无链表零树编码  designed by  陈军  2000-5-15
// 主程序
// 本程序在原LMZC程序的基础上吸收sazc的预扫描方式编制而成！
// 支持6144x6144，    2001－11－16

#include <stdio.h>
#include <memory.h>
#include <math.h>
#include <time.h> 
//#include "SpDecDll.h"
#include <xmmintrin.h>
#include <emmintrin.h>

#define WV_IMG_WIDTH 4016
#define	WV_IMG_HEIGHT	256
//Optimization [06.04.17]
#define WV_LEVEL     4
#define WV_SHIFT		2
#define WV_IMGSIZE0  WV_IMG_WIDTH*WV_IMG_HEIGHT

//Optimization [06.04.17]

#define ROUND(a) (((a)<0) ? (short) ((a)-0.5) : (short) ((a)+0.5))

int imgsize; //一棵生成树的大小，四级小波情况下为16*16
//int imgsize0;//整幅图像的大小
int imgsize1;
//     
short width, height;//一棵生成树的宽度、高度
short width0, height0, level;//原始图像的宽度、高度，小波分解级数
//                          
short *data,n;
short data0[WV_IMGSIZE0];
short coef[WV_IMGSIZE0];
short coef_int[WV_IMGSIZE0];
int tree_num;
//

int P_num[100],Sign_num[100],D_num[100],L_num[100];
int unitlen_x=1,unitlen_y=1,l=0,unitlen_x0,unitlen_y0;
int Fd_exist[256*4],Fl_exist[256*4],Fd_par_exist[256*4],Fl_par_exist[256*4],sn[256*4];
short thres;
int thres_ori;
//

short	d1[WV_IMG_WIDTH/2+2];
short	s1[WV_IMG_WIDTH/2+3];
short	dd[WV_IMG_WIDTH/2+4];
short	s[WV_IMG_WIDTH/2+3];
short	d0[WV_IMG_WIDTH/2];
short	s0[WV_IMG_WIDTH/2+1];


//__declspec(align(16)) short data0[WV_IMGSIZE0];
//__declspec(align(16)) short coef[WV_IMGSIZE0];
//__declspec(align(16)) short coef_int[WV_IMGSIZE0];
//__declspec(align(16)) short s1[4096+24] , s0[4096], d0[4096]; //大小取八列 与 一行二者中的最大值
//__declspec(align(16))__int64 const203[2],const217[2],const1817[2],const113[2],consts164[2],constd2048[2],consts064[2],constd12048[2],const_217[2] ,const_203[2];
//__declspec(align(16))__int64 const2048[2],const64[2];


#include "wv97wky.h"
#include "bitio.h"
#include "ispwky.h"






void megtree(int i,int j,int i0,int j0,int i1,int j1)
{
	int k0,l0,k1,l1,off;
	
	off=(i*unitlen_x0+j)*imgsize;

	for(k0=2*i0,k1=2*i1;k0<2*i0+2;k0++,k1++)
		for(l0=2*j0,l1=2*j1;l0<2*j0+2;l0++,l1++)
		{
			data0[k0*width0+l0]=coef[off+k1*width+l1];
			if(k0<height0/2 && l0<width0/2) 
				megtree(i,j,k0,l0,k1,l1);
		}
}
void megtee(int i,int j)
{
	int off;

	off = (i*unitlen_x0+j)*imgsize;

	data0[i*width0+j]                         = coef[off+0];
	data0[i*width0+j+unitlen_x0]              = coef[off+1];
	data0[(i+unitlen_y0)*width0+j]            = coef[off+width];
	data0[(i+unitlen_y0)*width0+j+unitlen_x0] = coef[off+width+1];

	megtree(i,j,i,j+unitlen_x0,0,1);
	megtree(i,j,i+unitlen_y0,j,1,0);
	megtree(i,j,i+unitlen_y0,j+unitlen_x0,1,1);
}
void megtee0()//  
{
	__m128i *pcoef,*pcoef1;
	__m128i *pdata0,*pdata01;
	int offp;//one tree has 32*8 point
	int half_img=WV_IMGSIZE0/16;
	int w = width0/8;
	int i,j,l;

	pcoef = (__m128i *)coef;
	pdata0 = (__m128i *)data0;

//--------------------LH0-----------------------------------
	for (j = 0; j < unitlen_y0; j++)
		for (i = 0; i < unitlen_x0; i ++)//(i,j)block
		{
			offp = (j*unitlen_x0+i)*32;
			for(l = 0; l < 8; l ++)//line in the block
			{
				pcoef1 = pcoef+offp+2*l+1;
				pdata01 = pdata0+(j*8+l)*w+w/2+i;
                *pdata01 = *pcoef1;
			}
		}

//--------------------HL0 ----------------------------------
	for (j = 0; j < unitlen_y0; j++)
		for (i = 0; i < unitlen_x0; i ++)//(i,j)block
		{
			offp = (j*unitlen_x0+i)*32;
			for(l = 0; l < 8; l ++)
			{
				pcoef1 = pcoef+offp+16+2*l;//16---half of a tree,
				pdata01 = pdata0+half_img+(j*8+l)*w+i;
				*pdata01 = *pcoef1;
			}
			
		}

		
//---------------------HH0---------------------------------------
	for (j = 0; j < unitlen_y0; j++)
		for (i = 0; i < unitlen_x0; i ++)//(i,j)block
		{
			offp = (j*unitlen_x0+i)*32;
			for(l = 0; l < 8; l ++)
			{
				pcoef1 = pcoef+offp+16+2*l+1;
				pdata01 = pdata0+half_img+(j*8+l)*w+w/2+i;
				*pdata01 = *pcoef1;
			}
		}

	
}




void megtee1()
{
	__int64 *pcoef,*pcoef1;
	__int64 *pdata0,*pdata01;
	int offp;
	int w = width0/4;
	int qua_img = WV_IMGSIZE0/16;
	int i,j,l;

	pcoef = (__int64 *)coef;
	pdata0 = (__int64 *)data0;

//------------------LH1---------------------------------------
	for (j = 0; j < unitlen_y0; j++)
		for (i = 0; i < unitlen_x0; i ++)//(i,j)block
		{
			offp = (j*unitlen_x0+i)*64;
			for(l = 0; l < 4; l ++)
			{
				pcoef1 = pcoef+offp+4*l+1;
				pdata01 = pdata0+(j*4+l)*w+w/4+i;
				*pdata01 = *pcoef1;
			}
		}

//------------------HL1--------------------------------------------
	for (j = 0; j < unitlen_y0; j++)
		for (i = 0; i < unitlen_x0; i ++)//(i,j)block
		{
			offp = (j*unitlen_x0+i)*64;
			for(l = 0; l < 4; l ++)
			{
				pcoef1 = pcoef+offp+16+4*l;
				pdata01 = pdata0+qua_img+(4*j+l)*w+i;
				*pdata01 = *pcoef1;
			}
		}
		
//------------------HH1------------------------------------------------
	for (j = 0; j < unitlen_y0; j++)
		for (i = 0; i < unitlen_x0; i ++)//(i,j)block
		{	
			offp = (j*unitlen_x0+i)*64;
			for(l = 0; l < 4; l ++)
			{
				pcoef1 = pcoef+offp+16+4*l+1;
				pdata01 = pdata0+qua_img+(4*j+l)*w+w/4+i;
				*pdata01 = *pcoef1;
			}
		}
		

}




void megtee2()
{
	__int32 *pcoef,*pcoef1;
	__int32 *pdata0,*pdata01;
	int offp;
	int w = width0/2;
	int hqua_img = WV_IMGSIZE0/16;
	int i,j,l;

	pcoef = (__int32 *)coef;
	pdata0 = (__int32 *)data0;

//----------------LH2----------------------------------------
	for (j = 0; j < unitlen_y0; j++)
		for (i = 0; i < unitlen_x0; i ++)//(i,j)block
		{
			offp = (j*unitlen_x0+i)*128;
			for(l = 0; l < 2; l ++)
			{
				pcoef1 = pcoef+offp+8*l+1;
				pdata01 = pdata0+(2*j+l)*w+w/8+i;
				*pdata01 = *pcoef1;
			}
		}
		
//----------------HL2------------------------------------------
	for (j = 0; j < unitlen_y0; j++)
		for (i = 0; i < unitlen_x0; i ++)//(i,j)block
		{
			offp = (j*unitlen_x0+i)*128;
			for(l = 0; l < 2; l ++)
			{
				pcoef1 = pcoef+offp+16+8*l;
				pdata01 = pdata0+hqua_img+(2*j+l)*w+i;
				*pdata01 = *pcoef1;
			}
		}

//----------------HH2---------------------------------------------
	for (j = 0; j < unitlen_y0; j++)
		for (i = 0; i < unitlen_x0; i ++)//(i,j)block
		{
			offp = (j*unitlen_x0+i)*128;
			for(l = 0; l < 2; l ++)
			{
				pcoef1 = pcoef+offp+16+8*l+1;
				pdata01 = pdata0+hqua_img+(2*j+l)*w+w/8+i;
				*pdata01 = *pcoef1;
			}
		}
	
}




void megtee3()
{
	__int16 *pcoef,*pcoef1;
	__int16 *pdata0,*pdata01;
	int w = width0;
	int offp;
	int i,j;

	pcoef = (__int16 *)coef;
	pdata0 = (__int16 *)data0;


//-----------LH3-----------------------------------------------
	for (j = 0; j < unitlen_y0; j++)
		for (i = 0; i < unitlen_x0; i ++)//(i,j)block
		{
			offp = (j*unitlen_x0+i)*256;
			pcoef1 = pcoef+offp+1;
			pdata01 = pdata0+j*w+unitlen_x0+i;
			*pdata01 = *pcoef1;
		}

//------------HL3-----------------------------------------------
	for (j = 0; j < unitlen_y0; j++)
		for (i = 0; i < unitlen_x0; i ++)//(i,j)block
		{
			offp = (j*unitlen_x0+i)*256;
			pcoef1 = pcoef+offp+16;
			pdata01 = pdata0+WV_IMGSIZE0/16+j*w+i;
			*pdata01 = *pcoef1;
		}


//-------------HH3-----------------------------------------------
	for (j = 0; j < unitlen_y0; j++)
		for (i = 0; i < unitlen_x0; i ++)//(i,j)block
		{
			offp = (j*unitlen_x0+i)*256;
			pcoef1 = pcoef+offp+16+1;
			pdata01 = pdata0+WV_IMGSIZE0/16+j*w+unitlen_x0+i;
			*pdata01 = *pcoef1;
		}
		
//-------------LL3-----------------------------------------------
	for (j = 0; j < unitlen_y0; j++)
		for (i = 0; i < unitlen_x0; i ++)//(i,j)block
		{
			offp = (j*unitlen_x0+i)*256;
			pcoef1 = pcoef+offp;
			pdata01 = pdata0+j*w+i;
			*pdata01 = *pcoef1;
		}


}


int fnSpDecDll(unsigned char *pchSpihtStr, 
		  unsigned char *pchReconImg,
		  int            nSpihtStrlen)
{
	int shift;
	int nDecLen;
	int i,k,l;
	int nTem=0;
// 	FILE *fpOut;

	width0   = WV_IMG_WIDTH; //整幅图像的宽度
	height0  = WV_IMG_HEIGHT; //整幅图像的高度
	level    = WV_LEVEL; //小波级数，默认是4;
	shift    = WV_SHIFT; //系数移位数，默认为2; 
	bitlen   = nSpihtStrlen;

	width    = 1<<level;   //一棵生成树的宽度
	height   = 1<<level;   //一棵生成树的高度
	
	unitlen_y0 = height0 >> level; //最低频LL子带的高度
	unitlen_x0 = width0  >> level; //最低频LL子带的宽度

//	imgsize0   = WV_IMG_WIDTH * height0; //整幅图像的大小
	imgsize    = width  * height;  //一棵生成树的大小，四级小波情况下为16*16
		
//	data0    = (short *)calloc(imgsize0,sizeof(short));//
//	coef     = (short *)calloc(imgsize0,sizeof(short));
//	coef_int = (short *)calloc(imgsize0,sizeof(short));	
	memset(coef,0,WV_IMGSIZE0*sizeof(short));

	//-----------------------Spiht decoder--------------------------------
//
//	if((fpOut=fopen("SpLib.dat","wb"))==NULL){
//		printf("Spdec.dat open error!\n");
//		return 0;
//	}
	bit = pchSpihtStr;

//fwrite(bit,sizeof(unsigned char),WV_IMGSIZE0,fpOut);
//fflush(fpOut);
	spiht_decode();
	nDecLen = bitlen;
//fwrite(coef,sizeof(unsigned char),WV_IMGSIZE0,fpOut);
//fflush(fpOut);
	for(k=0; k<unitlen_y0; k++)
		for(l=0; l<unitlen_x0; l++)
			megtee(k,l);
//	megtee0(); 	megtee1();
// 	megtee2();  megtee3();
//fwrite(data0,sizeof(unsigned char),WV_IMGSIZE0,fpOut);
//fflush(fpOut);
	//-----------------------Spiht decoder--------------------------------

	//-------------------invert wavelet transform-------------------------
//	printf("  (4) Doing invert wavelet transform\n");
	for(i=0; i<WV_IMGSIZE0; i++)
	{
		coef_int[i] = (data0[i]<<shift);
//		coef_int[i] = (data0[i]*(1<<shift));
	}

//	iwt(level, width0, height0, coef_int);
	iwt(coef_int);

	for (i=0; i<WV_IMGSIZE0; i++)
	{
		data0[i] = (coef_int[i]>>shift);
	}		
	
	for(i=0; i<WV_IMGSIZE0; i++)
	{
		nTem = data0[i];
		if(nTem>255)    nTem = 255;
		else if(nTem<0) nTem = 0;
		pchReconImg[i] = (unsigned char) nTem; 
	}
	//-------------------invert wavelet transform-------------------------

	//----------------------reconstruct image-----------------------------
//	printf("*** Spiht decoding finish!\n\n");
//	fprintf(g_fpDecRpt, "*** Spiht decoding finish!\n\n");

	return nDecLen;
}