sueipofi.c

1.Polarization analysis and filtering for three-component data 2.SUEIPOFI - EIgenimage (SVD) based

                
                /* construct and save eigenimages */
                for (j=0; j<3; j++) {
                    for (i=0; i<3; i++) {
                        for (it=jt; it<jt+iwl; it++) {
                            eigen[j][i][it] = \
                                udata[it-jt][j] * vdata[i][j] * sval[j];
                        }
                    }
                }

                /* calculate polarization attributes */
                for (j=0; j<2; j++) initpa[j][iw] = calc_svdrl(sval,j,pwr);
                initpa[2][iw] = calc_planarity(sval,pwr);

            } /* end loop time windows */

            /* interpolation of all polarization attributes */
            /* initially calculated values are asigned to window centers */
            for (i=0; i<3; i++) {
                if (iscubic) {
                    /* interpolation via cubic splines */
                    do_intcub( initpa[i], polar[i], nw, iwl, iwl/2, nt);
                    
                    /*... and clip values outside valid range */
                    for (it=0;it<nt; it++) {
                        if (polar[i][it]>1.0) polar[i][it] = 1.0;
                        if (polar[i][it]<0.0) polar[i][it] = 0.0;
                    }
                        
                }
                else {
                    /* linear interpolation */
                    do_intlin( initpa[i], polar[i], nw, iwl, iwl/2, nt);
                }
            }
            
            /* filter: polarization attribute weighted eigenimage stack */
            /* loop over components */

            for (i=0; i<3; i++) {
                /*loop over samples */
                for (it=0; it<(nw*iwl); it++) {
                    data[i][it] = eigen[0][i][it]*polar[0][it] + \
                                  eigen[1][i][it]*polar[1][it];
                    data[i][it] *= polar[2][it];
                }
            }

            /* write three-component data to stdout */
            fputdata3c(stdout, headerfp, data, nt);
            
            /* write polarization attributes to output files */
            if (rl1) fputdata(rl1fp, headerfp, polar[0], nt);
            if (rl2) fputdata(rl2fp, headerfp, polar[1], nt);
            if (pln) fputdata(plnfp, headerfp, polar[2], nt);            
            
        } /* end of processing three-component dataset */
        
    } while (gettr(&tr)); /* end loop over traces */
    
    if (verbose) {
        fprintf(stderr,"\n");
        if (icomp) warn("last %d trace(s) skipped", icomp);
    }
    
    /* close files for polarization attributes */
    if (rl1) efclose(rl1fp);
    if (rl2) efclose(rl2fp);
    if (pln) efclose(plnfp);
      
    return(CWP_Exit());
}

/**********************************************************************/
/* Functions used for polarization analysis                           */
/**********************************************************************/

/* rectilinearity, calculated from singular values                    */
/* (after Jurkevics, 1988, Franco and Musacchio, 2000)                */
/* w[3] = array containing singular values in descending order        */

float calc_svdrl(float *w, int axis, float pwr)
{
    float rl;   /* rectilinearity */

    if (w[axis]*w[axis]) {
        rl = 1.0 - (w[2]*w[2]) / (w[axis]*w[axis]);
    }
    else {
        rl = 0.0;
    }
        
    return pow(rl,pwr);
}

/* planarity, calculated from singular values (Jurkevics, 1988)       */
/* w[3] = array containing singular values in descending order        */
float calc_planarity(float *w, float pwr)
{
    float p;    /* planarity */
    
    if (w[0] + w[1]) {
        p = 1.0 - (2.0 * w[2]*w[2]) / (w[0]*w[0] + w[1]*w[1]);
    }
    else {
        p = 0.0;
    }
    
    return pow(p,pwr);
}

/**********************************************************************/
/* interpolation                                                      */
/**********************************************************************/

/* For more information on cubic spline interpolation methods used    */
/* here, see comments in "$CWPROOT/src/cwp/lib/cubicspline.c" and     */
/* "$CWPROOT/src/cwp/lib/intcub.c".                                   */

/* perform cubic spline interpolation of polarization attributes */
/* here: nin=nw, din=iwl, nout=nt, inx0=iwl/2                    */

void do_intcub(float *iny, float *outy, int nin, int din, int inx0, int nout)
{
    register int ix;
    float *inx, *outx;
    float (*yd)[4];
    
    /* allocate space */
    yd = (float(*)[4])ealloc1float(4*nin);
    inx = ealloc1float(nin);
    outx = ealloc1float(nout);

    /* initialize abscissa */   
    for (ix=0; ix<nin; ix++) {
        inx[ix] = (float) (inx0 + ix * din);
    }
    for (ix=0; ix<nout; ix++) {
        outx[ix] = (float) ix;
    }

    /* compute cubic spline coefficients */
    csplin(nin,inx,iny,yd);

    /* cubic spline interpolation */
    intcub(0,nin,inx,yd,nout,outx,outy);
    
    /* free memory */
    free(yd);
    free1float(inx);
    free1float(outx);

}

/* perform linear interpolation of polarization attributes       */
/* here: nin=nw, din=iwl, nout=nt, inx0=iwl/2                    */

void do_intlin(float *iny, float *outy, int nin, int din, int inx0, int nout)
{
    register int ix;
    float *inx, *outx;
    
    /* allocate space */
    inx = ealloc1float(nin);
    outx = ealloc1float(nout);

    /* initialize abscissa */   
    for (ix=0; ix<nin; ix++) {
        inx[ix] = (float) (inx0 + ix * din);
    }
    for (ix=0; ix<nout; ix++) {
        outx[ix] = (float) ix;
    }

    /* linear interpolation */
    intlin(nin,inx,iny,iny[0],iny[nin],nout,outx,outy);
    
    /* free memory */
    free1float(inx);
    free1float(outx);

}

/**********************************************************************/
/* utility functions                                                  */
/**********************************************************************/

/* write one-component data into file */

void fputdata(FILE *fileptr, FILE *headerptr, float *outdata, int nt)
{    
    efread(&tr, 1, HDRBYTES, headerptr);
    erewind(headerptr);
    
    memcpy((void *)tr.data, (const void *) outdata, nt*FSIZE);

    fputtr(fileptr, &tr);
}

/* write three-component data into file */

void fputdata3c(FILE *fileptr, FILE *headerptr, float **outdata3c, int nt)
{
    int i;
   
    for(i=0;i<3;i++) {
        efread(&tr, 1, HDRBYTES, headerptr);
        
        memcpy((void *)tr.data, (const void *) outdata3c[i], nt*FSIZE);

        fputtr(fileptr, &tr);
    }
    erewind(headerptr);
}
/* END OF FILE */