sopen_scp_read.c

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

	struct DATA_INFO textual;
	bool   AS_DECODE = 0; 

	decode.length_BdR0 = NULL; 	
	decode.samples_BdR0= NULL;
	decode.length_Res  = NULL;
	decode.samples_Res = NULL;
	decode.t_Huffman=NULL;
	decode.flag_Huffman=NULL;
	decode.data_lead=NULL;
	decode.data_protected=NULL;
	decode.data_subtraction=NULL;
	decode.length_BdR0=NULL;
	decode.samples_BdR0=NULL;
	decode.Median=NULL;
	decode.length_Res=NULL;
	decode.samples_Res=NULL;
	decode.Residual=NULL;
	decode.Reconstructed=NULL;

	//variables inizialization
	decode.flag_lead.number=0;
	decode.flag_lead.subtraction=0;
	decode.flag_lead.all_simultaneously=0;
	decode.flag_lead.number_simultaneously=0;

	decode.flag_BdR0.length=0;
	decode.flag_BdR0.fcM=0;
	decode.flag_BdR0.AVM=0;
	decode.flag_BdR0.STM=0;
	decode.flag_BdR0.number_samples=0;
	decode.flag_BdR0.encoding=0;

	decode.flag_Res.AVM=0;
	decode.flag_Res.STM=0;
	decode.flag_Res.number=0;
	decode.flag_Res.number_samples=0;
	decode.flag_Res.encoding=0;
	decode.flag_Res.bimodal=0;
	decode.flag_Res.decimation_factor=0;
#endif 
	
	ptr = hdr->AS.Header; 
	hdr->NRec = 0;

	sectionStart = 6;
	PtrCurSect = ptr+sectionStart;	

	/**** SECTION 0 ****/
	len = leu32p(PtrCurSect+4); 
	NSections = min((len-16)/10,_NUM_SECTION);

	section[0].ID	  = 0; 	
	section[0].length = len; 	
	section[0].index  = 6+16;
	for (int K=1; K<_NUM_SECTION; K++) {
		section[K].ID	  = K; 	
		section[K].length = 0; 	
		section[K].index  = 0;
	}

	for (int K=1; K<NSections; K++)	{
		// this is needed because fields are not always sorted
		curSect = leu32p(ptr+6+16+K*10); 
		if (curSect < NSections) {
			section[curSect].ID 	= curSect; 	
			section[curSect].length = leu32p(ptr+6+16+K*10+2); 	
			section[curSect].index  = leu32p(ptr+6+16+K*10+6)-1;
		}
	}

	for (int K=1; K<NSections; K++)	{
		curSect           = section[K].ID;
		len		  = section[K].length;
		sectionStart 	  = section[K].index;

	if (VERBOSE_LEVEL>8)
		fprintf(stdout,"SCP Section %i %i len=%i secStart=%i HeaderLength=%i\n",K,curSect,len,sectionStart,hdr->HeadLen);

	if (len==0) continue;	 /***** empty section *****/
		PtrCurSect = ptr+sectionStart;
		crc 	   = leu16p(PtrCurSect);
		uint16_t tmpcrc = CRCEvaluate((uint8_t*)(PtrCurSect+2),len-2); 
		if ((crc != 0xffff) && (crc != tmpcrc))
			fprintf(stderr,"Warning SOPEN(SCP-READ): faulty CRC in section %i: crc=%x, %x\n",curSect,crc,tmpcrc);
		if (curSect != leu16p(PtrCurSect+2))
			fprintf(stderr,"Warning SOPEN(SCP-READ): Current Section No does not match field in sections (%i %i)\n",curSect,leu16p(PtrCurSect+2)); 
		if (len != leu32p(PtrCurSect+4))
			fprintf(stderr,"Warning SOPEN(SCP-READ): length field in pointer section (%i) does not match length field in sections (%i %i)\n",K,len,leu32p(PtrCurSect+4)); 

		uint8_t versionSection  = *(ptr+sectionStart+8);
		uint8_t versionProtocol = *(ptr+sectionStart+9);

		curSectPos = 16;

		/**** SECTION 0 ****/
		if (curSect==0)  
		{
		}		

		/**** SECTION 1 ****/
		else if (curSect==1)  
		{
			struct tm t0,t1;
			t0.tm_year = 0;
			t0.tm_mon  = 0;
			t0.tm_mday = 0;
			t0.tm_hour = 0;
			t0.tm_min  = 0;
			t0.tm_sec  = 0;
			t0.tm_isdst= -1; // daylight savings time - unknown 
			hdr->T0    = 0;
			hdr->Patient.Birthday = 0;
			uint32_t len1;
 
 			while (*(PtrCurSect+curSectPos) < 255) {
				tag = *(PtrCurSect+curSectPos);
				len1 = leu16p(PtrCurSect+curSectPos+1);
				if (VERBOSE_LEVEL>8)
					fprintf(stdout,"SCP(r): Section 1 Tag %i Len %i\n",tag,len1);

				curSectPos += 3;
				if (curSectPos+len1 > len) {
					fprintf(stdout,"Warning SCP(read): section 1 corrupted (exceeds file length)\n");
			break;
				} 	 
				if (tag==0) {
					if (!hdr->FLAG.ANONYMOUS)
						strncpy(hdr->Patient.Name, (char*)(PtrCurSect+curSectPos),min(len1,MAX_LENGTH_NAME));
				}		
				else if (tag==1) {
//					hdr->Patient.FirstName = (char*)(PtrCurSect+curSectPos);
				}
				else if (tag==2) {
					if (len1>MAX_LENGTH_PID) {
						fprintf(stdout,"Warning SCP(read): length of Patient Id (section1 tag2) exceeds %i>%i\n",len1,MAX_LENGTH_PID); 
					}	
					strncpy(hdr->Patient.Id,(char*)(PtrCurSect+curSectPos),min(len1,MAX_LENGTH_PID));

					if (!strcmp(hdr->Patient.Id,"UNKNOWN")) 
						hdr->Patient.Id[0] = 0;
				}
				else if (tag==3) {
				}
				else if (tag==4) {
				}
				else if (tag==5) {
					t1.tm_year = leu16p(PtrCurSect+curSectPos)-1900;
					t1.tm_mon  = *(PtrCurSect+curSectPos+2)-1;
					t1.tm_mday = *(PtrCurSect+curSectPos+3);
					t1.tm_hour = 12; 
					t1.tm_min  =  0; 
					t1.tm_sec  =  0; 
					t1.tm_isdst= -1; // daylight saving time: unknown
//					t1.tm_gmtoff  =  0; 
					hdr->Patient.Birthday = tm_time2gdf_time(&t1);
				}
				else if (tag==6) {
					hdr->Patient.Height = leu16p(PtrCurSect+curSectPos);
				}
				else if (tag==7) {
					hdr->Patient.Weight = leu16p(PtrCurSect+curSectPos);
				}
				else if (tag==8) {
					hdr->Patient.Sex = *(PtrCurSect+curSectPos);
				}
				else if (tag==9) {
				}
				else if (tag==10) {
				}
				else if (tag==11) {
 					aECG->diastolicBloodPressure = leu16p(PtrCurSect+curSectPos);
				}
				else if (tag==12) {
					aECG->systolicBloodPressure  = leu16p(PtrCurSect+curSectPos);
				}
				else if (tag==13) {
					aECG->Diagnosis = (char*)(PtrCurSect+curSectPos);
				}
				else if (tag==14) {
					if (len1>80)
						fprintf(stderr,"Warning SCP(r): length of tag14 %i>40\n",len1);
					memcpy(hdr->ID.Manufacturer._field,(char*)PtrCurSect+curSectPos,min(len1,MAX_LENGTH_MANUF));
					hdr->ID.Manufacturer._field[min(len1,MAX_LENGTH_MANUF)] = 0;
					hdr->ID.Manufacturer.Model = hdr->ID.Manufacturer._field+8;  
					hdr->ID.Manufacturer.Version = hdr->ID.Manufacturer._field+36;  
					int tmp = strlen(hdr->ID.Manufacturer.Version)+1;
					hdr->ID.Manufacturer.SerialNumber = hdr->ID.Manufacturer.Version+tmp;
					tmp += strlen(hdr->ID.Manufacturer.Version+tmp)+1;	// skip SW ID
					tmp += strlen(hdr->ID.Manufacturer.Version+tmp)+1;	// skip SW
					tmp += strlen(hdr->ID.Manufacturer.Version+tmp)+1;	// skip SW
					hdr->ID.Manufacturer.Name = hdr->ID.Manufacturer.Version+tmp;
					
					/* might become obsolete */					
					//memcpy(hdr->aECG->Section1.tag14,PtrCurSect+curSectPos,40);
					//hdr->VERSION = *(PtrCurSect+curSectPos+14)/10.0;	// tag 14, byte 15
					aECG->Section1.Tag14.INST_NUMBER = leu16p(PtrCurSect+curSectPos);
					aECG->Section1.Tag14.DEPT_NUMBER = leu16p(PtrCurSect+curSectPos+2);
					aECG->Section1.Tag14.DEVICE_ID   = leu16p(PtrCurSect+curSectPos+4);
					aECG->Section1.Tag14.DEVICE_TYPE = *(PtrCurSect+curSectPos+ 6);
					aECG->Section1.Tag14.MANUF_CODE  = *(PtrCurSect+curSectPos+ 7);	// tag 14, byte 7 (MANUF_CODE has to be 255)
					aECG->Section1.Tag14.MOD_DESC    = (char*)(PtrCurSect+curSectPos+8); 
					aECG->Section1.Tag14.VERSION     = *(PtrCurSect+curSectPos+14);
					aECG->Section1.Tag14.PROT_COMP_LEVEL = *(PtrCurSect+curSectPos+15); 	// tag 14, byte 15 (PROT_COMP_LEVEL has to be 0xA0 => level II)
					aECG->Section1.Tag14.LANG_SUPP_CODE  = *(PtrCurSect+curSectPos+16);	// tag 14, byte 16 (LANG_SUPP_CODE has to be 0x00 => Ascii only, latin and 1-byte code)
					aECG->Section1.Tag14.ECG_CAP_DEV     = *(PtrCurSect+curSectPos+17);	// tag 14, byte 17 (ECG_CAP_DEV has to be 0xD0 => Acquire, (No Analysis), Print and Store)
					aECG->Section1.Tag14.MAINS_FREQ      = *(PtrCurSect+curSectPos+18);	// tag 14, byte 18 (MAINS_FREQ has to be 0: unspecified, 1: 50 Hz, 2: 60Hz)

					aECG->Section1.Tag14.ANAL_PROG_REV_NUM = (char*)(PtrCurSect+curSectPos+36);
					tmp = strlen((char*)(PtrCurSect+curSectPos+36));					
					aECG->Section1.Tag14.SERIAL_NUMBER_ACQ_DEV = (char*)(PtrCurSect+curSectPos+36+tmp+1);
					tmp += strlen((char*)(PtrCurSect+curSectPos+36+tmp+1));					
					aECG->Section1.Tag14.ACQ_DEV_SYS_SW_ID = (char*)(PtrCurSect+curSectPos+36+tmp+1);
					tmp += strlen((char*)(PtrCurSect+curSectPos+36+tmp+1));					
					aECG->Section1.Tag14.ACQ_DEV_SCP_SW = (char*)(PtrCurSect+curSectPos+36+tmp+1); 	// tag 14, byte 38 (SCP_IMPL_SW has to be "OpenECG XML-SCP 1.00")
					tmp += strlen((char*)(PtrCurSect+curSectPos+36+tmp+1));
					aECG->Section1.Tag14.ACQ_DEV_MANUF  = (char*)(PtrCurSect+curSectPos+36+tmp+1);	// tag 14, byte 38 (ACQ_DEV_MANUF has to be "Manufacturer")
				}
				else if (tag==15) {
					//memcpy(hdr->aECG->Section1.tag15,PtrCurSect+curSectPos,40);
					aECG->Section1.Tag15.VERSION     = *(PtrCurSect+curSectPos+14);
				}
				else if (tag==16) {
				}
				else if (tag==17) {
				}
				else if (tag==18) {
				}
				else if (tag==19) {
				}
				else if (tag==20) {
					aECG->ReferringPhysician = (char*)(PtrCurSect+curSectPos);
				}
				else if (tag==21) {
					aECG->MedicationDrugs = (char*)(PtrCurSect+curSectPos);
				}
				else if (tag==22) {
				}
				else if (tag==23) {
				}
				else if (tag==24) {
					aECG->EmergencyLevel = *(PtrCurSect+curSectPos);
				}
				else if (tag==25) {
					t0.tm_year = leu16p(PtrCurSect+curSectPos)-1900;
					t0.tm_mon  = (*(PtrCurSect+curSectPos+2)) - 1;
					t0.tm_mday = *(PtrCurSect+curSectPos+3);
				}
				else if (tag==26) {
					t0.tm_hour = *(PtrCurSect+curSectPos);
					t0.tm_min  = *(PtrCurSect+curSectPos+1);
					t0.tm_sec  = *(PtrCurSect+curSectPos+2);
				}
				else if (tag==27) {
					HighPass   = leu16p(PtrCurSect+curSectPos)/100.0;
				}
				else if (tag==28) {
					LowPass    = leu16p(PtrCurSect+curSectPos);
				}
				else if (tag==29) {
					uint8_t bitmap = *(PtrCurSect+curSectPos);
					if (bitmap==0)
						Notch = NaN;	// undefined 
					else if ((bitmap & 0x03)==0)
						Notch = -1;	// notch off
					else if (bitmap & 0x01)
						Notch = 60.0; 	// notch 60Hz
					else if (bitmap & 0x02)
						Notch = 50.0; 	// notch 50Hz
				}
				else if (tag==30) {
				}
 				else if (tag==31) {
				}
				else if (tag==32) {
				}
				else if (tag==33) {
				}
				else if (tag==34) {
					int16_t tzmin = lei16p(PtrCurSect+curSectPos);
					if (tzmin != 0x7fff)
						if (abs(tzmin)<=780) 
							t0.tm_min += tzmin;
						else 
							fprintf(stderr,"Warning SOPEN(SCP-READ): invalid time zone (Section 1, Tag34)\n");
					//fprintf(stdout,"SOPEN(SCP-READ): tzmin = %i %x \n",tzmin,tzmin);
				}
				else {
				}
				curSectPos += len1;
			}
//			t0.tm_gmtoff = 60*tzminutes;
			t0.tm_isdst = -1;
			hdr->T0     = tm_time2gdf_time(&t0);
		}

		/**** SECTION 2 ****/
		else if (curSect==2)  {
			aECG->FLAG.HUFFMAN = 1; 
			en1064.FLAG.HUFFMAN = 1; 

			NHT = leu16p(PtrCurSect+curSectPos);
			curSectPos += 2;
			if (NHT==19999) {
				en1064.FLAG.HUFFMAN = 1; 
				Huffman = (huffman_t*)malloc(sizeof(huffman_t));
				HTrees  = (htree_t**)malloc(sizeof(htree_t*));
				Huffman[0].NCT   = 19;
				Huffman[0].Table = DefaultTable;
				HTrees [0] = makeTree(Huffman[0]);
				k2 = 0; 
				if (VERBOSE_LEVEL==9)
					for (k1=0; k1<Huffman[k2].NCT; k1++) 
					fprintf(stdout,"%3i: %2i %2i %1i %3i %6u \n",k1,Huffman[k2].Table[k1].PrefixLength,Huffman[k2].Table[k1].CodeLength,Huffman[k2].Table[k1].TableModeSwitch,Huffman[k2].Table[k1].BaseValue,Huffman[k2].Table[k1].BaseCode); 
				if (!checkTree(HTrees[0])) // ### OPTIONAL, not needed ###
					fprintf(stderr,"Warning: invalid Huffman Tree\n");
			}
			else {
				en1064.FLAG.HUFFMAN = NHT; 
				Huffman = (huffman_t*)malloc(NHT*sizeof(huffman_t));
				for (k2=0; k2<NHT; k2++) {
					Huffman[k2].NCT   = leu16p(PtrCurSect+curSectPos);
					curSectPos += 2;
					Huffman[k2].Table = (typeof(Huffman[k2].Table))malloc(Huffman[k2].NCT * sizeof(*Huffman[k2].Table));
					HTrees      = (htree_t**)malloc(Huffman[k2].NCT*sizeof(htree_t*));
					for (k1=0; k1<Huffman[k2].NCT; k1++) {
						Huffman[k2].Table[k1].PrefixLength = *(PtrCurSect+curSectPos);  
						Huffman[k2].Table[k1].CodeLength = *(PtrCurSect+curSectPos+1);  
						Huffman[k2].Table[k1].TableModeSwitch = *(PtrCurSect+curSectPos+2);  
						Huffman[k2].Table[k1].BaseValue  = lei16p(PtrCurSect+curSectPos+3);  
						Huffman[k2].Table[k1].BaseCode   = leu32p(PtrCurSect+curSectPos+5);
						curSectPos += 9;
						if (VERBOSE_LEVEL==9)
							fprintf(stdout,"%3i %3i: %2i %2i %1i %3i %6u \n",k2,k1,Huffman[k2].Table[k1].PrefixLength,Huffman[k2].Table[k1].CodeLength,Huffman[k2].Table[k1].TableModeSwitch,Huffman[k2].Table[k1].BaseValue,Huffman[k2].Table[k1].BaseCode);
					}
					HTrees[k2] = makeTree(Huffman[k2]);
					if (!checkTree(HTrees[k2])) {
						B4C_ERRMSG = "Warning: invalid Huffman Tree\n";
						B4C_ERRNUM = B4C_DECOMPRESSION_FAILED; 
						// AS_DECODE = 2; // forced use of SCP-DECODE 
					}
				}
			}
		}

		/**** SECTION 3 ****/
		else if (curSect==3)  
		{
			hdr->NS = *(PtrCurSect+curSectPos);
			aECG->FLAG.REF_BEAT = (*(PtrCurSect+curSectPos+1) & 0x01);
			en1064.FLAG.REF_BEAT = (*(PtrCurSect+curSectPos+1) & 0x01);
			en1064.Section3.flags = *(PtrCurSect+curSectPos+1);
			if (en1064.FLAG.REF_BEAT && !section[4].length)
				fprintf(stderr,"Warning (SCP): Reference Beat but no Section 4\n");
			if (!(en1064.Section3.flags & 0x04) || ((en1064.Section3.flags>>3) != hdr->NS))
				fprintf(stderr,"Warning (SCP): channels are not simultaneously recorded! %x %i\n",en1064.Section3.flags,hdr->NS);

			curSectPos += 2;
			hdr->CHANNEL = (CHANNEL_TYPE *) calloc(hdr->NS, sizeof(CHANNEL_TYPE));
			en1064.Section3.lead = (typeof(en1064.Section3.lead))malloc(hdr->NS*sizeof(*en1064.Section3.lead));

			uint32_t startindex0; 
			startindex0 = leu32p(PtrCurSect+curSectPos);
			for (i = 0, hdr->SPR=1; i < hdr->NS; i++) {
				en1064.Section3.lead[i].start = leu32p(PtrCurSect+curSectPos);
				en1064.Section3.lead[i].end   = leu32p(PtrCurSect+curSectPos+4);

				hdr->CHANNEL[i].SPR 	= en1064.Section3.lead[i].end - en1064.Section3.lead[i].start + 1;
				hdr->SPR 		= lcm(hdr->SPR,hdr->CHANNEL[i].SPR);
				hdr->CHANNEL[i].LeadIdCode =  *(PtrCurSect+curSectPos+8);
				hdr->CHANNEL[i].Label[0]= 0;
				hdr->CHANNEL[i].LowPass = LowPass; 
				hdr->CHANNEL[i].HighPass= HighPass; 
				hdr->CHANNEL[i].Notch 	= Notch; 
				curSectPos += 9;

				if (en1064.Section3.lead[i].start != startindex0)
					fprintf(stderr,"Warning SCP(read): starting sample %i of #%i differ to %x in #1\n",en1064.Section3.lead[i].start,*(PtrCurSect+curSectPos+8),startindex0);
			}
		}
		/**** SECTION 4 ****/
		else if (curSect==4)  {