sudecomp.c

seismic software,very useful

char *sdoc =
"SUDECOMP - wavelet decompression of SU datasets in multiple parts\n"
"\n"
"sudecomp < file.comp [options]		report statistics of comp file\n"
"sudecomp < file.comp > file.su [options]\n"
"\n"
"options:\n"
"start=0 count=all end=end    decompress these blocks only\n"
"min=all max=all              decompress only values between this range\n"
"head=1                       assume SEGY file header is prepended to file\n"
"verbose=0                    report key value, size, and decompression speed\n"
"                             if > 1, report every verbose# chunks\n"
"datain=                      input compressed dataset name \n"
"                             (instead of standard input) \n"
"                             specify datain=file1 datain=file2 ... \n"
"                             datain=filen for multiple input data sets \n"
"skipfile=                    name of the file to indicate how many traces \n"
"                             to skip before output. The number of traces \n"
"                             to skip is writen starting at the scolumn on \n"
"                             the last line of the skipfile, as explained \n"
"                             below.  \n"
"scolumn=29                   starting colume position of the last line of \n"
"                             the skipfile to indicate number of traces \n"
"                             to skip \n"
"on2trace=0                   output total number of traces in input to \n"
"                             the key word tracl (first 4 bytes) in the \n"
"                             first output trace \n"
"                             0=no 1=yes \n"
"\n"
"Notes:\n"
"1) MSE is mean squared error of compression.\n"
"   It is non-zero if it was specified as a forward compression parameter.\n"
"\n";
#include <stdio.h>
#include <su.h>
#include <segy.h>
#include <cwp.h>
#include <hdr.h>
#include "sucomp.h"
#include "seilib.h"

struct comp_header ch;
segytrace tr;

main (int argc , char **argv) {
	int verbose=0, igroup=0, group1=0, group2=2000000000;
	int count=0, maxtrace=0, ntrace=0, ns=0, output=0;
	int hbytes=0, hflag=0, alltrace=0, ngroup=0;
	int min=-2000000000, max=2000000000, head=1;
	int hmin=2000000000, hmax=-2000000000;
	long maxbyte=0;
	long blen = 0;
	float *data=0, time0=0, alltime=0, timex=0, compress=0, mse=0;
	char *cbuf = (char *) NULL;
	FILE *outfp=stdout;
	float tmp;
	int itmp;

	FILE *infp; 
	char *datain;
	int ndata;
	int *nitrace, mtrace;
	char *skipfile;
	int ntraceskip, i, idata, idatain;
	int scolumn, on2trace, firsttrace, itrace, ntotal;

/* SU initializations */
	initargs (argc,argv);
	askdoc (1);

	if (filestat(1) != TTY) {
		output = 1;
		file2g(outfp);
	}

/* fetch parameters */
	getparint ("start",&group1);
	getparint ("end",&group2);
	if (getparint ("count",&count)) group2 = group1 + count - 1;
	getparint ("min",&min);
	getparint ("max",&max);
	getparint ("verbose",&verbose);
	getparint ("head",&head);
	if(!getparint("on2trace",&on2trace)) on2trace=0;

	ndata = 0;
	ndata = countparname("datain");
	idatain = 1;
	if(ndata==0) {
		ndata = 1;
		idatain = 0;
	}
	nitrace = (int*) malloc(ndata*sizeof(int));

/* scan input file */
	if (output) {fprintf (stderr,"preliminary pass through headers ...\n"); 
	fflush (stderr);}
	if (!output && verbose) fprintf (stderr,"Headers:\n");

	ngroup = 0;
	for(i=0;i<ndata;i++) {
		if(idatain==0) {
			infp = stdin;
		} else {
		  	getnparstring(i+1,"datain",&datain);
			fprintf (stderr," reading headers of %s\n",datain);
			infp = fopen(datain, "r");
		}
		file2g(infp);
		nitrace[i] = 0;
		if (head) fseeko(infp,3600,0);
		fread (&ch,1,sizeof(ch),infp);
		do {
			if (strcmp (ch.id,SUCOMPID)) {
				fprintf (stderr,"error at %d: %s\n",ngroup,ch.id);
				err("header error=%d loc=%d message=%s\n",
					ngroup,(int)fseeko(infp,0,1),ch.id);
			} else	{
				if (ch.ntrace < 0) err ("ch.ntrace group %d negative",ngroup);
				maxtrace = maxtrace > ch.ntrace ? maxtrace : ch.ntrace;
				alltrace += ch.ntrace;
				nitrace[i] += ch.ntrace;
				if (ch.nbyte < 0) err ("ch.nbyte group %d negative",ngroup);
				maxbyte = maxbyte > ch.nbyte ? maxbyte : ch.nbyte;
				hmin = hmin < vtoi(hdtype(ch.key),ch.value) ? hmin : vtoi(hdtype(ch.key),ch.value);
				hmax = hmax > vtoi(hdtype(ch.key),ch.value) ? hmax : vtoi(hdtype(ch.key),ch.value);
				compress += ch.compress;
				mse += ch.mse;
				if (!output && verbose && ngroup%verbose==0 
					&& ngroup>=group1 && ngroup<=group2 
					&& vtoi(hdtype(ch.key),ch.value) >= min 
					&& vtoi(hdtype(ch.key),ch.value) <= max) {
					fprintf (stderr, "%s=%d ntraces=%d bytes=%d compress=%g mse=%g\n", ch.key, vtoi(hdtype(ch.key),ch.value), ch.ntrace, ch.nbyte, ch.compress, ch.mse);
				}

			}
			ngroup += 1;
			if (fseeko(infp,ch.nbyte,1) < 0) break;

		} while( fread (&ch,1,sizeof(ch),infp) );

		if(idatain!=0) fprintf (stderr," There are %d traces in %s\n",
			nitrace[i],datain);
		if(idatain!=0) fclose(infp);
	}

	group2 = ngroup < group2 ? ngroup-1 : group2;
	ngroup = ngroup > 0 ? ngroup : 1;
	fprintf (stderr,"Compressed file summary:\n");
	fprintf (stderr,"key=%s min=%d max=%d\n",ch.key,hmin,hmax);
	fprintf (stderr,"ngroup=%d start=%d end=%d\n",ngroup,group1,group2);
	fprintf (stderr,"total_traces=%d largest_traces=%d largest_bytes=%d\n",alltrace, maxtrace, maxbyte);
	fprintf (stderr,
		"mean_compress=%g mean_mse=%g\n",compress/ngroup,mse/ngroup);

/* decompression */
	if (!output) exit (0);

/* find out the skip trace number */
	if (getparstring("skipfile",&skipfile)) {
		if(!getparint("scolumn",&scolumn)) scolumn=29;
		if(scolumn>80) err("scolumn must be less than 80");
		if( (infp = fopen(skipfile,"r"))!=NULL ) {
			cbuf = (char*) malloc(80*sizeof(char));
			for(i=0;i<2000000000;i++) {
				if(feof(infp) !=0) break;
				bzero(cbuf,80);
				fgets(cbuf,80,infp);
				sscanf(cbuf+scolumn-1,"%d",&ntraceskip);
			}
			fprintf(stderr," number of traces to skip = %d \n",ntraceskip);
			fclose(infp);
			free(cbuf);
		} else {
			ntraceskip = 0;
			fprintf(stderr," skipfile not present: set ntraceskip=0  \n");
			fprintf(stderr," number of traces to skip = %d \n",ntraceskip);
		}
	} else {
		ntraceskip = 0;
	}

	idata = ndata;
	itmp = 0;
	for(i=0;i<ndata;i++) {
		itmp += nitrace[i];
		if(itmp>ntraceskip) {
			idata = i;
			break;
		}
	}
	ntotal = alltrace;


/* allocate storage */
	data = (float*) malloc (maxbyte>3600?maxbyte:3600);
	alltime = 0;
	alltrace = 0;
	ngroup = 0;
	igroup = 0;
	mtrace = 0;
	for(i=0;i<idata;i++) mtrace += nitrace[i];
	firsttrace = 1;

	for(i=idata;i<ndata;i++) {
		if(idatain!=0) {
		  	getnparstring(i+1,"datain",&datain);
			infp = fopen(datain, "r");
		}
		file2g(infp);

		/* restore header */
		fseeko (infp,0,0);
		if (head) {
			fread (data,1,3600,infp);
			if(i==idata) fwrite (data,1,3600,outfp);
		}
		fread (&ch,1,sizeof(ch),infp);
/* decompress each chunk */
		do {
			if (   vtoi(hdtype(ch.key),ch.value) >= min 
			&& vtoi(hdtype(ch.key),ch.value) <= max
			&& mtrace+ch.ntrace>ntraceskip ) {
				if (verbose && igroup%verbose == 0) 
				fprintf (stderr, "#%d %s=%d ntraces=%d comp_bytes=%d compress=%g mse=%g", igroup, ch.key, vtoi(hdtype(ch.key),ch.value), ch.ntrace, ch.nbyte, ch.compress, ch.mse);

				time0 = cputime ();
				if (fread (data,1,ch.nbyte,infp) < ch.nbyte) break;
				cbuf = (char *) NULL;
				if(mtrace>=ntraceskip) {
					if(ch.ntrace>=18) {
						seismic_decompress ((char*)data,DATA_FORMAT,&ntrace,
							&ns,&hbytes,&hflag, ENDIAN,&blen,&cbuf);
						if(firsttrace==1 && on2trace==1 ) {
							firsttrace=0;
							bcopy(cbuf,&tr,240);
							tr.tracl = ntotal;
							bcopy(&tr,cbuf,240);
						}
						fwrite (cbuf,1,(ch.ns+60)*ch.ntrace*4,outfp);
					} else {
						if(firsttrace==1 && on2trace==1) {
							firsttrace=0;
							bcopy(data,&tr,240);
							tr.tracl = ntotal;
							bcopy(&tr,data,240);
						}
						fwrite (data,1,(ch.ns+60)*ch.ntrace*4,outfp);
					}
				} else {
					if(ch.ntrace>=18) {
						seismic_decompress ((char*)data,DATA_FORMAT,&ntrace,
							&ns,&hbytes,&hflag, ENDIAN,&blen,&cbuf);
						for(itrace=0;itrace<ch.ntrace;itrace++) {
							if(itrace+mtrace<ntraceskip) continue;
							bcopy(cbuf+itrace*(ch.ns+60)*4,&tr,(ch.ns+60)*4);
							if(firsttrace==1 && on2trace==1 ) {
								firsttrace=0;
								tr.tracl = ntotal;
							}
							fwrite (&tr,1,(ch.ns+60)*4,outfp);
						}
					} else {
						for(itrace=0;itrace<ch.ntrace;itrace++) {
							if(itrace+mtrace<ntraceskip) continue;
							bcopy(data+itrace*(ch.ns+60),&tr,(ch.ns+60)*4);
							if(firsttrace==1 && on2trace==1 ) {
								firsttrace=0;
								tr.tracl = ntotal;
							}
							fwrite (&tr,1,(ch.ns+60)*4,outfp);
						}
					}
				}
				if(cbuf != (char *) NULL ) free (cbuf);
				timex = cputime() - time0;
				if (verbose && igroup%verbose == 0) fprintf (stderr," sec=%g rate=%g\n",timex,(ch.ntrace*ch.ns*4)*.000001/timex);
				alltime += timex;
				alltrace += ch.ntrace;
				ngroup++;
			} else {
				fseeko (infp,ch.nbyte,1);
			}

			/* if unit in bounds */
			igroup += 1;
			if(igroup>group2) break;
			mtrace += ch.ntrace;

		} while(fread (&ch,1,sizeof(ch),infp));

		if(igroup>group2) break;
		fclose(infp);
	}


	tmp = alltrace;
	tmp = tmp * ch.ns;
	tmp = tmp / 250000.;
	itmp = tmp;
	if(alltime>0.) tmp = tmp /alltime;

	if (alltime > 0.) 
	fprintf (stderr,"%d %s's with %d traces and %d megabytes decompressed at %g megabytes per second\n", ngroup, ch.key, alltrace, itmp,tmp);
}