ezw.c

该程序把数字图像处理与小波变换结合起来

	printf("\n");
	//释放编码器占据的数据空间
	for (i=0; i< 4; i++){
		BasicCoderDealloc(DominantPassCoder[i]);
	}
	BasicCoderDealloc(SubordinateListCoder);
	free(DomList);
	free(SubList);
	//删除零树图占据的数据空间
	MapDealloc(map);
	SubbandSymbolDealloc(ztmap);
}

/*----------------------------------------------------------------------------*/	
/*----------------------------------------------------------------------------*/	
/* as described in "A Fast Technique For Identifying Zerotrees In 
 * The EZW Algorithm", ...
 */
void EZW2BuildZerotreeMap(WTRANSFORM *transform, SubbandSymbol *ztmap)
{
	int i, j;
	int n, temp;

	for (n=transform->nSubbands-1; n>=0; n--){
		if (n>=transform->nSubbands-3){
			//由最高频带开始
			for (j=0; j<transform->subbandVSize[n]; j++){
				for (i=0; i<transform->subbandHSize[n]; i++){

					temp = log2((int)fabs(transform->subbandPtr[n]
										[j*transform->subbandHSize[n]+i]));
					
					if (temp==0){
						ztmap->subbandPtr[n][j*ztmap->subbandHSize[n]+i] = 0;
					}
					else{
						ztmap->subbandPtr[n][j*ztmap->subbandHSize[n]+i] = 1<<temp;
					}
				}
			}
		}
		else if (n==0){
			//低频带部分
			for (j=0; j<transform->subbandVSize[0]; j++){
				for (i=0; i<transform->subbandHSize[0]; i++){

					temp = log2((int)fabs(transform->subbandPtr[0]
										[j*transform->subbandHSize[0]+i]));

					ztmap->subbandPtr[0][j*ztmap->subbandHSize[0]+i] 
						= 	  ztmap->subbandPtr[1][j*ztmap->subbandHSize[1] + i]
							| ztmap->subbandPtr[2][j*ztmap->subbandHSize[2] + i]
							| ztmap->subbandPtr[3][j*ztmap->subbandHSize[3] + i];

					if (temp!=0){
						ztmap->subbandPtr[0][j*ztmap->subbandHSize[0]+i] 
						= (1<<temp) 
							| ztmap->subbandPtr[0][j*ztmap->subbandHSize[0]+i];
					}
				}
			}
		}
		else{
			for (j=0; j<transform->subbandVSize[n]; j++){
				for (i=0; i<transform->subbandHSize[n]; i++){

					temp = log2((int)fabs(transform->subbandPtr[n]
										[j*transform->subbandHSize[n]+i]));
					
					ztmap->subbandPtr[n][j*ztmap->subbandHSize[n]+i] 
						= 	  ztmap->subbandPtr[n+3]
								[2*j*ztmap->subbandHSize[n+3] + 2*i]
							| ztmap->subbandPtr[n+3]
								[2*j*ztmap->subbandHSize[n+3] + 2*i+1]
							| ztmap->subbandPtr[n+3]
								[(2*j+1)*ztmap->subbandHSize[n+3] + 2*i]
							| ztmap->subbandPtr[n+3]
								[(2*j+1)*ztmap->subbandHSize[n+3] + 2*i+1];
					
					if (temp!=0){
						ztmap->subbandPtr[n][j*ztmap->subbandHSize[n]+i] 
						= (1<<temp)
							| ztmap->subbandPtr[n][j*ztmap->subbandHSize[n]+i];
							
					}
				}
			}
			
		}
	}
}

/*----------------------------------------------------------------------------*/	
/*----------------------------------------------------------------------------*/	
void EncodeDominantPassMortonScan(DLNode *DomList,SLNode *SubList,
											 int *SubListSize, WTRANSFORM *transform, 
											 MAP *map, int threshold, 
											 ArithEncoder *encoder, 
											 BasicCoder **DominantPassCoder, 
											 SubbandSymbol *ztmap, int ByteBudget)
{
	int i, j, m, n;
	DLNode node;
	char NodeStatus;
	Boolean HighestSubband=FALSE;
	int IZCount;
	BasicCoder *Coder;
	Boolean ParentSignificant;
	Boolean PreviousCoeffSignificant;
	Boolean CurrentCoeffSignificant = FALSE;
	int DomListSize=0, DomListIndex;

	/* initialize LL_0 subband */
	for (j=0; j<transform->subbandVSize[0]; j+=2){
		for (i=0; i<transform->subbandHSize[0]; i+=2){
			/* begin quad_coeff */
			for (m=0; m<4; m++){
				node.x = i + ScanOrderX[m];
				node.y = j + ScanOrderY[m];
				node.scale = 0;
				node.orientation = 0;
				
				NodeStatus = MapGetNodeSymbol(map, node.scale, node.orientation, 
															node.x, node.y);
				
				/* output code */
				if (!NodeStatus){
					/* not a significant coeff in the previous dominant pass */
					PreviousCoeffSignificant = CurrentCoeffSignificant;
					
					CurrentCoeffSignificant = FALSE;
					EncodeNode(SubList, SubListSize, transform, &node, ztmap, threshold);						
					
					if (node.code == POS || node.code == NEG){
						MapSetNodeSymbol(map, node.scale, node.orientation, 
							node.x, node.y, 1);
						CurrentCoeffSignificant = TRUE;
					}
					
					if (ROUND(ArithEncoderNBitsOutput(encoder)/8.0) > ByteBudget-1){
						EndEncoding = TRUE;
						return;
					}
					
					/* perform context switching */
					ParentSignificant = BASEBAND_PARENT_SIGNIFICANT;
					
					if (ParentSignificant && PreviousCoeffSignificant){
						Coder = DominantPassCoder[0];
					}
					else if (ParentSignificant){
						Coder = DominantPassCoder[1];
					}
					else if (PreviousCoeffSignificant){
						Coder = DominantPassCoder[2];
					}
					else{
						Coder = DominantPassCoder[3];
					}
					
					BasicCoderEncode(Coder, encoder, node.code, TRUE);
					
				}
				else{
				/* no code since the node is already significant in previous pass */
					CurrentCoeffSignificant = TRUE;
				}
				
				if (node.code!=ZTR || NodeStatus){
					/* insert the offsprings - do the LH_1 first */
					/* LH_0 subband in the main list */
					DomList[DomListSize].x = node.x;
					DomList[DomListSize].y = node.y;
					DomList[DomListSize].scale = 1;
					DomList[DomListSize].orientation = 0;
					DomListSize++;
				}
			}
			/* end quad_coeff */
		}
	}
	
	/* add the HL_1 and HH_1 etries into the list */
	m = DomListSize;

	/* HL_1 */
	for (i=0; i<m; i++){
		DomList[DomListSize].x = DomList[i].x;
		DomList[DomListSize].y = DomList[i].y;
		DomList[DomListSize].scale = 1;
		DomList[DomListSize].orientation = 1;
		DomListSize++;
	}
	
	/* HH_1 */
	for (i=0; i<m; i++){
		DomList[DomListSize].x = DomList[i].x;
		DomList[DomListSize].y = DomList[i].y;
		DomList[DomListSize].scale = 1;
		DomList[DomListSize].orientation = 2;
		DomListSize++;
	}

	/* start the scanning */
	DomListIndex = 0;
	do{
		if (DomListIndex < DomListSize){
			node = DomList[DomListIndex];
		
			NodeStatus = MapGetNodeSymbol(map, node.scale, node.orientation, 
									node.x, node.y);
			
			if (!NodeStatus){
				PreviousCoeffSignificant = CurrentCoeffSignificant;
				CurrentCoeffSignificant = FALSE;
				EncodeNode(SubList, SubListSize, transform, &node, ztmap, threshold);

				if (node.code == POS || node.code == NEG){
					MapSetNodeSymbol(map, node.scale, node.orientation, 
											node.x, node.y, 1);
					CurrentCoeffSignificant = TRUE;
				}

				if (ROUND(ArithEncoderNBitsOutput(encoder)/8.0) > ByteBudget-1){
					EndEncoding = TRUE;
					return;
				}

				
				if ((node.scale == transform->nsteps) && !HighestSubband){
					/* prevent entering the code again */
					HighestSubband = TRUE;
					/* switch to 3 symbols */
					for (i=0; i<4;i++){
						IZCount = ContextGetProb(DominantPassCoder[i]->context, IZ);
						ContextPutValue(DominantPassCoder[i]->context, -IZCount, IZ);
						ContextPutValue(DominantPassCoder[i]->context, IZCount, ZTR);
					}
				}
				

				/* check if parent is significant by looking at the map */
				/* special treatment for scale 1 subbands */
				if (node.scale==1){
					if (MapGetNodeSymbol(map, 0, 0, node.x, node.y)){
						ParentSignificant = TRUE; 
					}
					else{
						ParentSignificant = FALSE;
					}
				}
				else{
					if (MapGetNodeSymbol(map, node.scale - 1, node.orientation, 
												node.x>>1, node.y>>1)){
						ParentSignificant = TRUE;
					}
					else{
						ParentSignificant = FALSE;
					}
				}
				
				/* perform context switching */
				if (ParentSignificant && PreviousCoeffSignificant){
					Coder = DominantPassCoder[0];
				}
				else if (ParentSignificant){
					Coder = DominantPassCoder[1];
				}
				else if (PreviousCoeffSignificant){
					Coder = DominantPassCoder[2];
				}
				else{
					Coder = DominantPassCoder[3];
				}

				BasicCoderEncode(Coder, encoder, node.code, TRUE);
				
			}
			else{
				CurrentCoeffSignificant = TRUE;
			}

			
			/* check if need to insert children */
			if ((node.code != ZTR || NodeStatus) && node.scale<transform->nsteps){
				i = node.x<<1;
				j = node.y<<1;
				node.scale = node.scale+1;
				
				/* same orientation */
				for (m=0; m<4; m++){
					DomList[DomListSize].x = i+ScanOrderX[m];
					DomList[DomListSize].y = j+ScanOrderY[m];
					DomList[DomListSize].scale = node.scale;
					DomList[DomListSize].orientation = node.orientation;
					DomListSize++;
				}
			}

			DomListIndex++;
		}
		else{
			break;
		}		
	}while(1);

	return;
}


/*----------------------------------------------------------------------------*/	
/*----------------------------------------------------------------------------*/	
void EncodeNode(SLNode *SubList, int *SubListSize, WTRANSFORM *transform, 
					 DLNode *node, SubbandSymbol *ztmap, int threshold)
{	
	double value;
	int idx;
	static i=0;
	
	if (node->scale==0){
		idx=0;
	}
	else{
		idx = 3*node->scale - 2 + node->orientation;
	}

	value = transform->subbandPtr[idx]
				[node->y*transform->subbandHSize[idx] + node->x];

	/* NOTE:
	 * We need to use integer conparison. Without rounding to integer,
	 * there may be some coefficients that are less than threshold by 
	 * 0.5, but are not being coded as significant. 
	 * Previous version of code ignore this and obtain less number of
	 * significant coefficients.
	 */
	if (ROUND(fabs(value))>=threshold){
		/*
		if (ROUND(fabs(value))==threshold){
			printf("%d\n", i++);
		}
		*/
		if (value>0.0){
			node->code = POS;
		}
		else{
			node->code = NEG;
		}
		
		/* put node into subordinate list */
		SubList[*SubListSize].x = node->x;
		SubList[*SubListSize].y = node->y;
		SubList[*SubListSize].scale = node->scale;
		SubList[*SubListSize].orientation = node->orientation;
		SubList[*SubListSize].qvalue = (int)(ROUND(fabs(value)));
		SubList[*SubListSize].rvalue = 3*(threshold>>1);
		(*SubListSize)++;

		/* set magnitude to 0 */
		transform->subbandPtr[idx]
				[node->y*transform->subbandHSize[idx] + node->x] = 0.0;
	}
	else{
		if ((ztmap->subbandPtr[idx][node->y*ztmap->subbandHSize[idx]+node->x]
				& threshold)==0){
			node->code = ZTR;
		}
		else{
			node->code = IZ;
		}
	}
}

/*----------------------------------------------------------------------------*/	
/*----------------------------------------------------------------------------*/	
void EncodeSubordinatePass(SLNode *SubList, int SubListSize, int threshold, 
							ArithEncoder *encoder, BasicCoder *SubordinateListCoder, 
							int ByteBudget)
{
	int i;

	if (threshold>0){
		for (i=0; i<SubListSize; i++){
			
			if (ROUND(ArithEncoderNBitsOutput(encoder)/8.0) > ByteBudget-1){
				EndEncoding = TRUE;
				break;
			}
			
			if ((SubList[i].qvalue & threshold) != 0){ 
				/* output '1' */
				BasicCoderEncode(SubordinateListCoder, encoder, 1, TRUE);
				SubList[i].rvalue += (threshold>>1);
			}
			else{
				/* output '0' */
				BasicCoderEncode(SubordinateListCoder, encoder, 0, TRUE);
				SubList[i].rvalue -= (threshold>>1);
			}
		}
	}
}

/*----------------------------------------------------------------------------*/	
/*----------------------------------------------------------------------------*/	
void EZW2Error(char *fmt, ...)
{
	va_list argptr;
	
	va_start( argptr, fmt );
	fprintf(stderr, "EZW2Error: " );
	vprintf( fmt, argptr );
	va_end( argptr );
	exit( -1 );
}

/*----------------------------------------------------------------------------*/	
/*----------------------------------------------------------------------------*/	
void EZW2Warning(char *fmt, ...)
{
	va_list argptr;
	
	va_start( argptr, fmt );
	fprintf( stderr, "EZW2Warning: " );
	vprintf( fmt, argptr );
	va_end( argptr );
}