sopen_scp_read.c

c++编写的用于生物信号处理的软件库

/*

    $Id: sopen_scp_read.c,v 1.55 2008/03/14 08:30:43 schloegl Exp $
    Copyright (C) 2005,2006,2007 Alois Schloegl <[e-mail]>

    This file is part of the "BioSig for C/C++" repository 
    (biosig4c++) at http://biosig.sf.net/ 


    This program is free software; you can redistribute it and/or
    modify it under the terms of the GNU General Public License
    as published by the Free Software Foundation; either version 3
    of the License, or (at your option) any later version.


 */


// #define WITHOUT_SCP_DECODE    // use SCP-DECODE if needed, Bimodal, reference beat  

/*
	the experimental version needs a few more thinks: 
	- Bimodal and RefBeat decoding do not work yet

	- validation and testing  
*/


#include <stdio.h>             // system includes
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <float.h>
#include <fcntl.h>
#include <sys/stat.h>
#include <time.h>
#include <ctype.h>

#include "../biosig-dev.h"

#include "structures.h"
static const U_int_S _NUM_SECTION=12U;	//consider first 11 sections of SCP
static bool add_filter=true;             // additional filtering gives better shape, but use with care

#ifndef WITHOUT_SCP_DECODE
int scp_decode(HDRTYPE*, pointer_section*, DATA_DECODE&, DATA_RECORD&, DATA_INFO&, bool&);
#endif

// Huffman Tables         	
uint16_t NHT; 	/* number of Huffman tables */
struct table_t {
	uint8_t PrefixLength;
       	uint8_t CodeLength;
       	uint8_t TableModeSwitch;
       	int16_t BaseValue;
       	uint32_t BaseCode;
}; 
struct huffman_t {
      	uint16_t NCT; 	/* number of Code structures in Table #1 */
       	table_t *Table; 
} *Huffman;

struct htree_t {
	htree_t* child0;
	htree_t* child1;
	uint16_t idxTable;
} **HTrees;

table_t DefaultTable[19] = {
	{ 1,  1, 1, 0, 0 },
	{ 3,  3, 1, 1, 1 },
	{ 3,  3, 1,-1, 5 },
	{ 4,  4, 1, 2, 3 },
	{ 4,  4, 1,-2, 11},
	{ 5,  5, 1, 3, 7 },
	{ 5,  5, 1,-3, 23},
	{ 6,  6, 1, 4, 15},
	{ 6,  6, 1,-4, 47},
	{ 7,  7, 1, 5, 31},
	{ 7,  7, 1,-5, 95},
	{ 8,  8, 1, 6, 63},
	{ 8,  8, 1,-6, 191},
	{ 9,  9, 1, 7, 127},
	{ 9,  9, 1,-7, 383},
	{10, 10, 1, 8, 255},
	{10, 10, 1,-8, 767},
	{18, 10, 1, 0, 511},
	{26, 10, 1, 0, 1023}
};

/*
	This structure defines the fields used for "Annotated ECG" 
 */
struct en1064_t {
	char*		test;		/* test field for annotated ECG */
	
	float		diastolicBloodPressure;				
	float		systolicBloodPressure;	
	char*		MedicationDrugs;
	char*		ReferringPhysician;
	char*		LatestConfirmingPhysician;
	char*		Diagnosis;
	uint8_t		EmergencyLevel; /* 0: routine 1-10: increased emergency level */

	float		HeartRate;	
	float		P_wave[2]; 	/* start and end  */
	float		QRS_wave[2]; 	/* start and end  */
	float		T_wave[2]; 	/* start and end  */
	float		P_QRS_T_axes[3];

	/***** SCP only fields *****/
	struct {	
		uint8_t	HUFFMAN;	
		uint8_t	REF_BEAT;
		uint8_t	DIFF;// OBSOLETE
		uint8_t	BIMODAL;// OBSOLETE
	} FLAG;
        struct {
		//uint8_t tag14[41],tag15[41];
	        struct {
			uint16_t INST_NUMBER;		/* tag 14, byte 1-2  */
			uint16_t DEPT_NUMBER;		/* tag 14, byte 3-4  */
			uint16_t DEVICE_ID;		/* tag 14, byte 5-6  */
			uint8_t  DEVICE_TYPE;		/* tag 14, byte 7: 0: Cart, 1: System (or Host)  */
			uint8_t MANUF_CODE;		/* tag 14, byte 8 (MANUF_CODE has to be 255) */
			char*   MOD_DESC;		/* tag 14, byte 9 (MOD_DESC has to be "Cart1") */
			uint8_t VERSION;		/* tag 14, byte 15 (VERSION has to be 20) */
			uint8_t PROT_COMP_LEVEL;	/* tag 14, byte 16 (PROT_COMP_LEVEL has to be 0xA0 => level II) */
			uint8_t LANG_SUPP_CODE;		/* tag 14, byte 17 (LANG_SUPP_CODE has to be 0x00 => Ascii only, latin and 1-byte code) */
			uint8_t ECG_CAP_DEV;		/* tag 14, byte 18 (ECG_CAP_DEV has to be 0xD0 => Acquire, (No Analysis), Print and Store) */
			uint8_t MAINS_FREQ;		/* tag 14, byte 19 (MAINS_FREQ has to be 0: unspecified, 1: 50 Hz, 2: 60Hz) */
			char 	reserved[22]; 		/* char[35-19] reserved; */			
			char* 	ANAL_PROG_REV_NUM;
			char* 	SERIAL_NUMBER_ACQ_DEV;
			char* 	ACQ_DEV_SYS_SW_ID;
			char* 	ACQ_DEV_SCP_SW; 	/* tag 14, byte 38 (SCP_IMPL_SW has to be "OpenECG XML-SCP 1.00") */
			char* 	ACQ_DEV_MANUF;		/* tag 14, byte 38 (ACQ_DEV_MANUF has to be "Manufacturer") */
        	} Tag14, Tag15; 
        } Section1;
        struct {
        } Section2;
        struct {
        	uint8_t NS, flags;
        	struct { 
        		uint32_t start;
        		uint32_t end;
//        		uint8_t  id;
        	} *lead;
        } Section3;
        struct {
        	uint16_t len_ms, fiducial_sample, N;
       		uint32_t SPR;
        	struct { 
        		uint16_t btyp;
        		uint32_t SB;
        		uint32_t fcM;
        		uint32_t SE;
        		uint32_t QB;
        		uint32_t QE;
        	} *beat;
        } Section4;
        struct {
        	size_t   StartPtr;
        	size_t	 Length;
        	uint16_t AVM, dT_us;
        	uint8_t  DIFF; //diff: see FLAG 
        	uint16_t *inlen;
        	int32_t  *datablock;
        } Section5;
        struct {
        	size_t   StartPtr;
        	size_t	 Length;
        	uint16_t AVM, dT_us;
        	uint8_t  DIFF, BIMODAL; //diff, bimodal: see FLAG 
        	uint16_t *inlen;
        	int32_t  *datablock;
        } Section6;
} en1064;

/* new node in Huffman tree */
htree_t* newNode() {
	htree_t* T  = (htree_t*) malloc(sizeof(htree_t));
	T->child0   = NULL;
	T->child1   = NULL;
	T->idxTable = 0; 
	return(T);
}

/* check Huffman tree */
int checkTree(htree_t *T) {
	int v,v1,v2,v3; 
	
	v1 = (T->child0 == NULL) && (T->child0 == NULL) && (T->idxTable > 0);
	v2 = (T->idxTable == 0) && (T->child0 != NULL) && checkTree(T->child0);
	v3 = (T->idxTable == 0) && (T->child1 != NULL) && checkTree(T->child1);
	v = v1 || v2 || v3;
	if (!v) fprintf(stderr,"Warning: Invalid Node in Huffman Tree: %i %p %p\n",T->idxTable,T->child0,T->child1);
	return(v);
}

/* convert Huffman Table into a Huffman tree */
htree_t* makeTree(huffman_t HT) {
	uint16_t k1,k2;
	htree_t* T = newNode();
	htree_t* node;
	for (k1=0; k1<HT.NCT; k1++) {
		node = T; 
		uint32_t bc = HT.Table[k1].BaseCode;
		for (k2=0; k2<HT.Table[k1].CodeLength; k2++, bc>>=1) {
			if (bc & 0x00000001) {
				if (node->child1==NULL) node->child1 = newNode();
				node = node->child1;
			}
			else {
				if (node->child0==NULL) node->child0 = newNode();
				node = node->child0;
			}
		}
		node->idxTable = k1+1;
	}
	return(T); 
}

/* get rid of Huffman tree */
void freeTree(htree_t* T) {
	if (T->child0 != NULL) freeTree(T->child0);
	if (T->child1 != NULL) freeTree(T->child1);
	free(T);
}

int DecodeHuffman(htree_t *HTrees[], huffman_t *HuffmanTables, uint8_t* indata, size_t inlen, int32_t* outdata, size_t outlen) {
	uint16_t ActualTable = 0; 
	htree_t *node; 
	size_t k1, k2, i;
	uint32_t acc; 
	int8_t dlen,k3,r;

	k1=0, k2=0;
	node = HTrees[ActualTable];
	r = 0; i = 0;
	while ((k1 < inlen*8) && (k2 < outlen)) {
		r = k1 % 8; 
		i = k1 / 8;

		if (!node->idxTable) {
			if (indata[i] & (1<<(7-r))) {
				if (node->child1 != NULL)
					node = node->child1;
				else {
					B4C_ERRMSG = "Empty node in Huffman table! Do not know what to do !\n";
					B4C_ERRNUM = B4C_DECOMPRESSION_FAILED;					
					return(-1);
				}	
			}	
			else {
				if (node->child0 != NULL)
					node = node->child0;
				else {
					B4C_ERRMSG = "Empty node in Huffman table! Do not know what to do !\n";
					B4C_ERRNUM = B4C_DECOMPRESSION_FAILED;					
					return(-1);
				}	
			}	
			++k1;
		}

		r = k1 % 8; 
		i = k1 / 8;

		if (node->idxTable) {
			// leaf of tree reached
			table_t TableEntry = HuffmanTables[ActualTable].Table[node->idxTable - 1]; 
			dlen = TableEntry.PrefixLength - TableEntry.CodeLength;
			if (!TableEntry.TableModeSwitch)
				// switch Huffman Code 
				ActualTable = TableEntry.BaseValue;
			else if (dlen) {
				// no compression
				acc = 0;  //(uint32_t)(indata[i]%(1<<r)); 
				for (k3=0; k3*8-r < dlen; k3++)
					acc = (acc<<8)+(uint32_t)indata[i+k3];
				
				outdata[k2] = (acc >> (k3*8 - r - dlen)) & ((1L << dlen) - 1L) ; 
				if (outdata[k2] >= (1<<dlen-1))
					outdata[k2] -= 1<<dlen;
				k1 += dlen; 
				++k2;
			}
			else {
				// lookup Huffman Table 
				outdata[k2++] = TableEntry.BaseValue; 
			}
			// reset node to root
			node = HTrees[ActualTable];
		}
	}
	return(0);
};

void deallocEN1064(en1064_t en1064) {	
	/* free allocated memory */ 
	if (en1064.FLAG.HUFFMAN) {
		for (size_t k1=0; k1<en1064.FLAG.HUFFMAN; k1++) {
		 	if (NHT!=19999) free(Huffman[k1].Table);
		 	freeTree(HTrees[k1]); 
		} 
		free(Huffman);
		free(HTrees); 
	}	

	if (en1064.Section3.lead != NULL) 	free(en1064.Section3.lead);
	if (en1064.Section4.beat != NULL) 	free(en1064.Section4.beat);
	if (en1064.Section5.inlen != NULL) 	free(en1064.Section5.inlen);
	if (en1064.Section5.datablock != NULL) 	free(en1064.Section5.datablock);
	if (en1064.Section6.inlen != NULL) 	free(en1064.Section6.inlen);
//	if (en1064.Section6.datablock != NULL) 	free(en1064.Section6.datablock);
}

int sopen_SCP_read(HDRTYPE* hdr) {	
/*
	this function is a stub or placeholder and need to be defined in order to be useful. 
	It will be called by the function SOPEN in "biosig.c"

	Input: 
		char* Header	// contains the file content
		
	Output: 
		HDRTYPE *hdr	// defines the HDR structure accoring to "biosig.h"
*/	

	uint8_t*	ptr; 	// pointer to memory mapping of the file layout
	uint8_t*	PtrCurSect;	// point to current section 
	uint8_t*	Ptr2datablock=NULL; 	// pointer to data block 
	int32_t* 	data;		// point to rawdata
	uint16_t	curSect=0; 	// current section
	uint32_t 	len; 
	uint16_t 	crc; 
	uint32_t	i,k1,k2; 
	size_t		curSectPos;
	size_t 		sectionStart; 
	int 		NSections = 12;
	uint8_t		tag;
	float 		HighPass=0, LowPass=1.0/0.0, Notch=-1; 	// filter settings
	uint16_t	Cal5=0, Cal6=0, Cal0;	// scaling coefficients 
	uint16_t 	dT_us = 1000; 	// sampling interval in microseconds

	/* 
	   Try direct conversion SCP->HDR to internal data structure
		+ whole data is loaded once, then no further File I/O is needed. 
		- currently Huffman and Bimodal compression is not supported. 
	*/	

	aECG_TYPE* aECG;
	if (hdr->aECG == NULL) {
		hdr->aECG = malloc(sizeof(aECG_TYPE));
		aECG = (aECG_TYPE*)hdr->aECG;
		aECG->diastolicBloodPressure=0.0;				 
		aECG->systolicBloodPressure=0.0;
		aECG->MedicationDrugs = "\0";
		aECG->ReferringPhysician="\0";
		aECG->LatestConfirmingPhysician="\0";
		aECG->Diagnosis="\0";
		aECG->EmergencyLevel=0;
		hdr->ID.Technician = "nobody";
	}
	else 
		aECG = (aECG_TYPE*)hdr->aECG;

	
	aECG->Section1.Tag14.VERSION = 0; // acquiring.protocol_revision_number 
	aECG->Section1.Tag15.VERSION = 0; // analyzing.protocol_revision_number
	aECG->FLAG.HUFFMAN  = 0; 
	aECG->FLAG.DIFF     = 0; 
	aECG->FLAG.REF_BEAT = 0; 
	aECG->FLAG.BIMODAL  = 0;

	en1064.Section4.len_ms	 = 0;
	
#ifndef WITHOUT_SCP_DECODE
	struct pointer_section section[_NUM_SECTION];
	struct DATA_DECODE decode;
	struct DATA_RECORD record;