1 /*
2
3 Copyright (C) 2000,2005,2007-2020 Alois Schloegl <alois.schloegl@gmail.com>
4 This file is part of the "BioSig for C/C++" repository
5 (biosig4c++) at http://biosig.sf.net/
6
7
8 This program is free software; you can redistribute it and/or
9 modify it under the terms of the GNU General Public License
10 as published by the Free Software Foundation; either version 3
11 of the License, or (at your option) any later version.
12
13 This program is distributed in the hope that it will be useful,
14 but WITHOUT ANY WARRANTY; without even the implied warranty of
15 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 GNU General Public License for more details.
17
18 You should have received a copy of the GNU General Public License
19 along with this program. If not, see <http://www.gnu.org/licenses/>.
20
21 */
22
23 #include <assert.h>
24 #include <ctype.h>
25 #include <math.h>
26 #include <stdio.h>
27 #include <stdlib.h>
28 #include <string.h>
29 #include <strings.h>
30 #include <time.h>
31 #include "biosig.h"
32
33 #define min(a,b) (((a) < (b)) ? (a) : (b))
34 #define max(a,b) (((a) > (b)) ? (a) : (b))
35
36 #ifdef __cplusplus
37 extern "C" {
38 #endif
39 int savelink(const char* filename);
40 #ifdef __cplusplus
41 }
42 #endif
43
44 #ifdef WITH_PDP
45 void sopen_pdp_read(HDRTYPE *hdr);
46 #endif
47
48 extern const uint16_t GDFTYP_BITS[];
49
50
compare_uint16(const void * a,const void * b)51 int compare_uint16 (const void* a, const void* b) {
52 return ( (int)*(uint16_t*)a - (int)*(uint16_t*)b );
53 }
54
main(int argc,char ** argv)55 int main(int argc, char **argv) {
56
57 HDRTYPE *hdr;
58 size_t count, k1, ne=0;
59 char *source, *dest, *tmpstr;
60 biosig_options_type biosig_options;
61 biosig_options.free_text_event_limiter = ""; // default is no delimiter
62
63 enum FileFormat SOURCE_TYPE; // type of file format
64 struct {
65 enum FileFormat TYPE;
66 float VERSION;
67 } TARGET;
68 TARGET.TYPE=GDF;
69 TARGET.VERSION=-1.0;
70
71 FILE* fid=stdout;
72 int COMPRESSION_LEVEL=0;
73 int status;
74 uint16_t k;
75 uint16_t chansel = 0;
76 int TARGETSEGMENT=1; // select segment in multi-segment file format EEG1100 (Nihon Kohden)
77 int VERBOSE = 0;
78 char FLAG_ANON = 1;
79 char FLAG_CSV = 0;
80 char FLAG_JSON = 0;
81 char FLAG_DYGRAPH = 0;
82 char *argsweep = NULL;
83 double t1=0.0, t2=1.0/0.0;
84 uint16_t *CHANLIST = NULL;
85 uint16_t LenChanList = 0; // number of entries in CHANLIST
86
87 #ifdef CHOLMOD_H
88 char *rrFile = NULL;
89 int refarg = 0;
90 #endif
91
92 source = NULL;
93 dest = NULL;
94
95 for (k=1; k<argc; k++) {
96 if (!strcmp(argv[k],"-v") || !strcmp(argv[k],"--version") ) {
97 fprintf(stderr,"biosig2gdf (BioSig4C++) v%i.%i.%i\n", BIOSIG_VERSION_MAJOR, BIOSIG_VERSION_MINOR, BIOSIG_PATCHLEVEL);
98 fprintf(stderr,"Copyright (C) 2006-2020 by Alois Schloegl and others\n");
99 fprintf(stderr,"This file is part of BioSig http://biosig.sf.net - the free and\n");
100 fprintf(stderr,"open source software library for biomedical signal processing.\n\n");
101 fprintf(stderr,"BioSig is free software; you can redistribute it and/or modify\n");
102 fprintf(stderr,"it under the terms of the GNU General Public License as published by\n");
103 fprintf(stderr,"the Free Software Foundation; either version 3 of the License, or\n");
104 fprintf(stderr,"(at your option) any later version.\n\n");
105 }
106
107 else if (!strcmp(argv[k],"-h") || !strcmp(argv[k],"--help") ) {
108 fprintf(stderr,"\nusage: biosig2gdf [OPTIONS] SOURCE\n");
109 fprintf(stderr,"\nusage: biosig2gdf [OPTIONS] SOURCE DEST\n");
110 fprintf(stderr," SOURCE is the source file \n");
111 fprintf(stderr," SOURCE can be also network file bscs://<hostname>/ID e.g. bscs://129.27.3.99/9aebcc5b4eef1024 \n");
112 fprintf(stderr," DEST is the destination file \n");
113 fprintf(stderr,"\n Supported OPTIONS are:\n");
114 fprintf(stderr," -v, --version\n\tprints version information\n");
115 fprintf(stderr," -h, --help \n\tprints this information\n");
116 #ifdef CHOLMOD_H
117 fprintf(stderr," -r, --ref=MM \n\trereference data with matrix file MM. \n\tMM must be a 'MatrixMarket matrix coordinate real general' file.\n");
118 #endif
119 fprintf(stderr, " -a, --anon=yes (default)\n"
120 "\tanonymized data processing - personalize data (name, and birthday) is not processed but ignored.\n"
121 "\tThe patient can be still identified with the unique patient identifier (and an external database).\n"
122 "\tThis is for many cases sufficient (e.g. for research etc.). This mode can be turn off with\n"
123 " -n, --anon=no\n"
124 "\tThis will process personal information like name and birthday. One might want to use this mode\n"
125 "\twhen converting personalized patient data and no unique patient identifier is available.\n"
126 "\tIt's recommended to pseudonize the data, or to use the patient identifier instead of patient name and birthday.\n"
127 );
128 fprintf(stderr, " -c <chanlist> (default: 0)\n"
129 "\tThis flag is in preparation and not working yet.\n"
130 "\t<chanlist> is a comma-separated list of channel numbers.\n"
131 "\tIf <chanlist> contains any number smaller than 1 or larger than the number \n"
132 "\t of available channels, all channels are loaded.\n"
133 );
134 fprintf(stderr," --free-text-event-limiter=\";\"\n\tfree text of events limited to first occurrence of \";\" (only EDF+/BDF+ format)\n");
135 fprintf(stderr," -s=#\tselect target segment # (in the multisegment file format EEG1100)\n");
136 fprintf(stderr," -SWEEP=ne,ng,ns\n\tsweep selection of HEKA/PM files\n\tne,ng, and ns select the number of experiment, the number of group, and the sweep number, resp.\n");
137 fprintf(stderr," -VERBOSE=#, verbosity level #\n\t0=silent [default], 9=debugging\n");
138 fprintf(stderr," --chan=CHAN\n\tselect channel CHAN (0: all channels, 1: first channel, etc.)\n");
139 fprintf(stderr,"\n\n");
140 return(0);
141 }
142
143 else if (!strncmp(argv[k],"-VERBOSE",2)) {
144 VERBOSE = argv[k][strlen(argv[k])-1]-48;
145 #ifndef NDEBUG
146 // then VERBOSE_LEVEL is not a constant but a variable
147 VERBOSE_LEVEL = VERBOSE;
148 #endif
149 }
150 else if (!strncasecmp(argv[k],"-SWEEP=",7))
151 argsweep = argv[k]+6;
152
153 else if (!strcasecmp(argv[k],"-a") || !strcasecmp(argv[k],"--anon") )
154 FLAG_ANON = 1;
155
156 else if (!strcasecmp(argv[k],"-n") || !strcasecmp(argv[k],"--anon=no") )
157 FLAG_ANON = 0;
158
159 else if (!strcasecmp(argv[k],"-c") ) {
160 char *CS = strdup(argv[++k]);
161 char *CS0 = CS;
162 CHANLIST = realloc(CHANLIST,sizeof(uint16_t)*strlen(CS));
163 LenChanList = 0;
164 while (isdigit(*CS)) {
165 CHANLIST[LenChanList++] = strtol(CS,&CS,10);
166 if (*CS==0) break;
167 if (*CS != ',') {
168 LenChanList = 0;
169 break;
170 }
171 CS++;
172 }
173 CHANLIST[LenChanList]=0xffff;
174 free(CS0);
175 qsort(CHANLIST, LenChanList, sizeof(uint16_t), compare_uint16);
176 }
177 else if (!strncmp(argv[k],"--free-text-event-limiter=",26))
178 biosig_options.free_text_event_limiter = strstr(argv[k],"=") + 1;
179
180 #ifdef CHOLMOD_H
181 else if ( !strncmp(argv[k],"-r=",3) || !strncmp(argv[k],"--ref=",6) ) {
182 // re-referencing matrix
183 refarg = k;
184 }
185 #endif
186
187 else if (!strncmp(argv[k],"-s=",3)) {
188 TARGETSEGMENT = atoi(argv[k]+3);
189 }
190
191 else if (argv[k][0]=='[' && argv[k][strlen(argv[k])-1]==']' && (tmpstr=strchr(argv[k],',')) ) {
192 t1 = strtod(argv[k]+1,NULL);
193 t2 = strtod(tmpstr+1,NULL);
194 if (VERBOSE_LEVEL>7) fprintf(stderr,"[%f,%f]\n",t1,t2);
195 }
196
197 else {
198 break;
199 }
200
201 if (VERBOSE_LEVEL>7)
202 fprintf(stderr,"%s (line %i): biosig2gdf: arg%i = <%s>\n",__FILE__,__LINE__, k, argv[k]);
203
204 } // end of for (k=1; k<argc; k++)
205
206 switch (argc - k) {
207 case 0:
208 fprintf(stderr,"biosig2gdf: missing file argument\n");
209 fprintf(stderr,"usage: biosig2gdf [options] SOURCE\n");
210 fprintf(stderr," for more details see also biosig2gdf --help \n");
211 exit(-1);
212 case 1:
213 source = argv[k];
214 break;
215 case 2:
216 source = argv[k];
217 dest = argv[k+1];
218 fid = fopen(dest,"w");
219 break;
220 default:
221 fprintf(stderr,"biosig2gdf: extra arguments %d-%d (%s, etc) ignored\n",argc,k,argv[k+1]);
222 }
223
224 if (VERBOSE_LEVEL<0) VERBOSE=1; // default
225 if (VERBOSE_LEVEL>7) fprintf(stderr,"%s (line %i): BIOSIG2GDF %s started \n",__FILE__,__LINE__, source);
226
227 tzset();
228 hdr = constructHDR(0,0);
229 // hdr->FLAG.OVERFLOWDETECTION = FlagOverflowDetection;
230 hdr->FLAG.UCAL = 1;
231 hdr->FLAG.TARGETSEGMENT = TARGETSEGMENT;
232 hdr->FLAG.ANONYMOUS = FLAG_ANON;
233
234 if (argsweep) {
235 k = 0;
236 do {
237 hdr->AS.SegSel[k++] = strtod(argsweep+1, &argsweep);
238 } while (argsweep[0]==',' && (k < 5) );
239 }
240 hdr = sopen_extended(source, "r", hdr, &biosig_options);
241
242 #ifdef WITH_PDP
243 if (hdr->AS.B4C_ERRNUM) {
244 biosigERROR(hdr, 0, NULL); // reset error
245 sopen_pdp_read(hdr);
246 }
247 #endif
248 // HEKA2ITX hack
249 if (TARGET.TYPE==ITX) {
250 if (hdr->TYPE==HEKA) {
251 // hack: HEKA->ITX conversion is already done in SOPEN
252 dest = NULL;
253 }
254 else {
255 fprintf(stderr,"error: only HEKA->ITX is supported - source file is not HEKA file");
256 biosigERROR(hdr, B4C_UNSPECIFIC_ERROR, "error: only HEKA->ITX is supported - source file is not HEKA file");
257 }
258 }
259
260 if ((status=serror2(hdr))) {
261 destructHDR(hdr);
262 exit(status);
263 }
264
265 t1 *= hdr->SampleRate / hdr->SPR;
266 t2 *= hdr->SampleRate / hdr->SPR;
267 if (isnan(t1)) t1 = 0.0;
268 if (t2+t1 > hdr->NRec) t2 = hdr->NRec - t1;
269
270 if ( ( t1 - floor (t1) ) || ( t2 - floor(t2) ) ) {
271 fprintf(stderr,"ERROR BIOSIG2GDF: cutting from parts of blocks not supported; t1 (%f) and t2 (%f) must be a multiple of block duration %f\n", t1,t2,hdr->SPR / hdr->SampleRate);
272 biosigERROR(hdr, B4C_UNSPECIFIC_ERROR, "blocks must not be split");
273 }
274
275 if ((status=serror2(hdr))) {
276 destructHDR(hdr);
277 exit(status);
278 }
279 sort_eventtable(hdr);
280
281 /*******************************************
282 Channel selection
283 *******************************************/
284 uint16_t NS2=0; // number of channels (w/o hidden channels)
285
286 for (k=0; k < hdr->NS; k++) {
287 NS2 += (hdr->CHANNEL[k].OnOff > 0);
288 }
289 for (k = 0; k < LenChanList; k++) {
290 // check whether all arguments specify a valid channel number
291 if ((CHANLIST[k] < 1) || (NS2 < CHANLIST[k])) {
292 LenChanList = 0; // in case of an invalid channel, all channels are selected
293 break;
294 }
295 }
296
297 if (LenChanList > 0) {
298 fprintf(stderr, "argument -c <chanlist> is currently not supported, the argument is ignored ");
299 LenChanList=0;
300 }
301 if (LenChanList > 0) {
302 // all entries in argument -c specify valid channels
303 int chan=0;
304 int k1 = 0;
305 for (k=0; k<hdr->NS; k++)
306 if (hdr->CHANNEL[k].OnOff > 0) {
307 chan++; // count
308 while (CHANLIST[k1] < chan) k1++; // skip double entries in CHANLIST
309
310 if (CHANLIST[k1]==chan) {
311 hdr->CHANNEL[k].OnOff = 1;
312 k1++;
313 fprintf(stderr,"-- %d %d %d \n",k,k1,chan);
314 }
315 else if (CHANLIST[k1] > chan)
316 hdr->CHANNEL[k].OnOff = 0;
317 }
318 }
319
320 for (k=0; k<hdr->NS; k++) {
321 if ( (hdr->CHANNEL[k].OnOff > 0) && hdr->CHANNEL[k].SPR ) {
322 if ((hdr->SPR/hdr->CHANNEL[k].SPR)*hdr->CHANNEL[k].SPR != hdr->SPR)
323 fprintf(stderr,"Warning: channel %i might be decimated!\n",k+1);
324 };
325 }
326 if (CHANLIST) {free(CHANLIST); CHANLIST=NULL;}
327
328 hdr->FLAG.OVERFLOWDETECTION = 0;
329 hdr->FLAG.UCAL = 1;
330 hdr->FLAG.ROW_BASED_CHANNELS = 1;
331
332 if (VERBOSE_LEVEL>7)
333 fprintf(stderr,"%s (line %i): SREAD [%f,%f].\n",__FILE__,__LINE__,t1,t2);
334
335 if (hdr->NRec <= 0) {
336 // in case number of samples is not known
337 count = sread(NULL, t1, (size_t)-1, hdr);
338 t2 = count;
339 }
340 else {
341 if (t2+t1 > hdr->NRec) t2 = hdr->NRec - t1;
342 count = sread(NULL, t1, t2, hdr);
343 }
344
345 biosig_data_type* data = hdr->data.block;
346
347 if ((status=serror2(hdr))) {
348 destructHDR(hdr);
349 exit(status);
350 };
351
352 if (hdr->FILE.OPEN) {
353 sclose(hdr);
354 free(hdr->AS.Header);
355 hdr->AS.Header = NULL;
356 if (VERBOSE_LEVEL>7) fprintf(stderr,"%s (line %i): file closed\n",__FILE__,__LINE__);
357 }
358
359 SOURCE_TYPE = hdr->TYPE;
360
361 hdr->TYPE = GDF;
362 hdr->VERSION = 3.0;
363
364 hdr->FILE.COMPRESSION = COMPRESSION_LEVEL;
365
366 /*******************************************
367 convert to unique sampling rate and unique data type
368 *******************************************/
369
370 uint16_t gdftyp=0;
371 if (1) {
372 uint16_t nbits=0;
373 int flag_fnbits=0;
374 int flag_snbits=0;
375 int flag_unbits=0;
376 int flag_signed=0;
377
378 for (k=0; k < hdr->NS; k++) {
379 if (hdr->CHANNEL[k].OnOff && hdr->CHANNEL[k].SPR) {
380 hdr->CHANNEL[k].SPR = hdr->SPR;
381
382 uint16_t gt = hdr->CHANNEL[k].GDFTYP;
383 if ((15<gt) && (gt<20))
384 if (flag_fnbits < GDFTYP_BITS[gt]) flag_fnbits=GDFTYP_BITS[gt];
385 if (((gt<15) && !(gt%2)) || (511<gt))
386 if (flag_unbits < GDFTYP_BITS[gt]) flag_unbits=GDFTYP_BITS[gt];
387 if (((gt<15) && (gt%2)) || ((255<gt) && (gt<511)))
388 if (flag_snbits < GDFTYP_BITS[gt]) flag_snbits=GDFTYP_BITS[gt];
389 }
390 }
391
392 if (flag_snbits>flag_unbits) { nbits=flag_snbits; flag_signed=1;}
393 else if (flag_snbits<flag_unbits) { nbits=flag_unbits; flag_signed=0;}
394 else {nbits = flag_unbits+1; flag_signed=1;}
395
396 if (flag_fnbits==32) gdftyp=16;
397 else if (flag_fnbits==64) gdftyp=17;
398 else if (flag_fnbits==128) gdftyp=18;
399 else if (nbits>64) gdftyp=8-flag_signed;
400 else if (nbits>32) gdftyp=8-flag_signed;
401 else if (nbits>16) gdftyp=6-flag_signed;
402 else if (nbits>8) gdftyp=4-flag_signed;
403 else gdftyp=2-flag_signed;
404
405 hdr->SPR = 1;
406 hdr->NRec = hdr->data.size[1];
407 typeof(hdr->NS) NS=0;
408 for (k=0; k<hdr->NS; k++) {
409 if (hdr->CHANNEL[k].OnOff) {
410 NS++;
411 hdr->CHANNEL[k].GDFTYP = gdftyp;
412 hdr->CHANNEL[k].SPR = hdr->SPR;
413 }
414 }
415 hdr->AS.bpb = NS*GDFTYP_BITS[gdftyp];
416 }
417
418 /*********************************
419 Write data
420 *********************************/
421
422 hdr->FLAG.ANONYMOUS = FLAG_ANON;
423
424 /*
425 keep header data from previous file, in might contain optional data
426 (GDF Header3, EventDescription, hdr->SCP.SectionX, etc. )
427 and might still be referenced and needed.
428 */
429 void *tmpmem = hdr->AS.Header;
430 hdr->AS.Header = NULL;
431 size_t filesize=0;
432
433 // open_gdfwrite(hdr,stdout);
434 // encode header
435 hdr->TYPE=GDF;
436 hdr->VERSION=3.0;
437 struct2gdfbin(hdr);
438
439 // write header into buffer
440 filesize += fwrite (hdr->AS.Header, 1, hdr->HeadLen, fid);
441
442 count = 0;
443 char BLK[0x10000];
444 unsigned spb = (0x10000 / hdr->AS.bpb) * hdr->AS.bpb / GDFTYP_BITS[gdftyp];
445 while (count < hdr->data.size[0]*hdr->data.size[1]) {
446 if (count+spb > hdr->data.size[0]*hdr->data.size[1])
447 spb = hdr->data.size[0] * hdr->data.size[1] - count;
448
449 switch (gdftyp) {
450 case 1: {
451 int8_t *data = (void*)&BLK;
452 for (int k=0; k < spb; k++) data[k]=(int8_t)hdr->data.block[k+count];
453 break;
454 }
455 case 2: {
456 uint8_t *data = (void*)&BLK;
457 for (int k=0; k < spb; k++) data[k]=(int8_t)hdr->data.block[k+count];
458 break;
459 }
460 case 3: {
461 int16_t *data = (void*)&BLK;
462 for (int k=0; k < spb; k++) lei16a((int16_t)hdr->data.block[k+count], data+k);
463 break;
464 }
465 case 4: {
466 uint16_t *data = (void*)&BLK;
467 for (int k=0; k < spb; k++) leu16a((uint16_t)hdr->data.block[k+count], data+k);
468 break;
469 }
470 case 5: {
471 int32_t *data = (void*)&BLK;
472 for (int k=0; k < spb; k++) lei32a((int32_t)hdr->data.block[k+count], data+k);
473 break;
474 }
475 case 6: {
476 uint32_t *data = (void*)&BLK;
477 for (int k=0; k < spb; k++) leu32a((uint32_t)hdr->data.block[k+count], data+k);
478 break;
479 }
480 case 7: {
481 int64_t *data = (void*)&BLK;
482 for (int k=0; k < spb; k++) lei64a((int64_t)hdr->data.block[k+count], data+k);
483 break;
484 }
485 case 8: {
486 uint64_t *data = (void*)&BLK;
487 for (int k=0; k < spb; k++) leu64a((uint64_t)hdr->data.block[k+count], data+k);
488 break;
489 }
490 case 16: {
491 float *data = (void*)&BLK;
492 for (int k=0; k < spb; k++) lef32a((float)hdr->data.block[k+count], data+k);
493 break;
494 }
495 case 17: {
496 double *data = (void*)&BLK;
497 for (int k=0; k < spb; k++) lef64a((double)hdr->data.block[k+count], data+k);
498 break;
499 }
500 default: {
501 }
502 }
503 count += fwrite(BLK, GDFTYP_BITS[gdftyp]/8, spb, fid);
504 }
505
506 filesize += count*GDFTYP_BITS[gdftyp]/8;
507
508 // write event table into buffer
509 size_t len3 = hdrEVT2rawEVT(hdr);
510 filesize += fwrite(hdr->AS.rawEventData, 1, len3, fid);
511 if (fid != stdout) fclose(fid);
512
513 status = serror2(hdr);
514 destructHDR(hdr);
515 exit(status);
516 }
517
518