suxarb.c

seismic software,very useful

      }
      fprintf( stderr ,"\n" );
   }

   /*-----------------------*/
   /* open the output files */
   /*-----------------------*/

   for( i=0; i<nFiles; i++ ){
      if( !(outfp[i]=fopen( outfile[i] ,"w" )) ){
         sprintf(  errString 
                   ,"fopen(3c) failed %s %d" ,__FILE__ ,__LINE__ );
         err( errString );
      }

   }

   /*---------------------*/
   /* read in first trace */
   /*---------------------*/

   if (!fgettr(stdin ,&tr) ){
      err("can't get first trace");
   }

   xtype  = hdtype(xkey);
   xindex = getindex(xkey);

   ytype  = hdtype(ykey);
   yindex = getindex(ykey);

   iltype  = hdtype(ilkey);
   ilindex = getindex(ilkey);

   xltype  = hdtype(xlkey);
   xlindex = getindex(xlkey);
   
   /*----------------------------*/
   /* process all the input data */
   /*----------------------------*/

   do {
   
      gethval(&tr, xindex, &xval);
      x0 = vtof(xtype, xval);

      gethval(&tr, yindex, &yval);
      y0 = vtof(ytype, yval);

      if( pass == 0 ){

         x0 *= dx;
         y0 *= dy;

         /*---------------------------------*/
         /* polyline transform & extraction */
         /*---------------------------------*/

         nl = PolylineTransform( x0 
                           ,y0 
                           ,xpoly[0] 
                           ,ypoly[0] 
                           ,spoly[0]
                           ,npoly[0]
                           ,il
                           ,xl
                          );

         for( j=0; j<nl; j++ ){

            if( fabs(xl[j]) <= dw ){

               changeval( iltype ,&ilval ,(int) il[j] );
               puthval( &tr ,ilindex ,&ilval );

               changeval( xltype ,&xlval ,(int) xl[j] );
               puthval( &tr ,xlindex ,&xlval );
            
               fputtr( outfp[0] ,&tr );
            }

         }

      }else{

         /*--------------------*/
         /* polygon extraction */
         /*--------------------*/

         if( pass == 1 ){

            for( i=0; i<nPolygons; i++ ){

               inside = InPolygon( x0 
                                  ,y0 
                                  ,xpoly[i] 
                                  ,ypoly[i] 
                                  ,npoly[i] 
                                 );

                if( inside ){
                  fputtr( outfp[i] ,&tr );
  
                }

            }

         }else if( pass == 2 ){

            for( i=0; i<nPolygons; i++ ){

               inside = InPolygon( x0 
                                  ,y0 
                                  ,xpoly[i] 
                                  ,ypoly[i] 
                                  ,npoly[i] 
                                 );

               if( inside ){
                  goto next;
               }

            }

            fputtr( outfp[0] ,&tr );

         }

      }

      next:;
         
   } while (fgettr(stdin ,&tr));


   return( 0 );
}

/*--------------------------------------------------------------------*\
   InPolygon() determines if a point is within a polygon denoted by
   a polyline with identical first and last points.

   The basic algorithm is from "Computational Geometry in C", by
   Joseph O'Rourke but has been completely rewritten to make it 
   suitable for industrial use.

   The function returns integer values indicating the location of
   the input point relative to the bounding polygon:

   0 - point is exterior to the polygon
   1 - point is a vertex of the polygon
   2 - point lies on a line segment
   3 - point is strictly interior to the polygon

   Reginald H. Beardsley                            [e-mail]
\*--------------------------------------------------------------------*/

#include <stdio.h>

int InPolygon( float X0 ,float Y0 ,float* X ,float* Y ,int n ){

   int i;            /* vertex          */
   int j;            /* adjacent vertex */

   int R = 0;        /* number of right crossings */
   int L = 0;        /* number of left crossings  */

   float x;

   for( i=0; i<n; i++ ){

      if( X[i] == X0 && Y[i] == Y0 ){
         return( 1 );

      }

      j = (i + n - 1) % n;

      if(   ((Y[i] > Y0) && (Y[j] <= Y0)) 
          ||((Y[j] > Y0) && (Y[i] <= Y0)) ){
         
         x = ((X[i]-X0) * (Y[j]-Y0) - (X[j]-X0) * (Y[i]-Y0))
            / ((Y[j]-Y0) - (Y[i]-Y0));

         if( x > 0 ){
            R++;
         }
      }

      if(   ((Y[i] > Y0) && (Y[j] <= Y0)) 
          ||((Y[j] > Y0) && (Y[i] <= Y0)) ){
         
         x = ((X[i]-X0) * (Y[j]-Y0) - (X[j]-X0) * (Y[i]-Y0))
            / ((Y[j]-Y0) - (Y[i]-Y0));

         if( x < 0 ){
            L++;
         }
      }

   }

   if( (R % 2) != (L % 2) ){
      return( 2 ); 
   }

   if( (R % 2) == 1 ){
      return( 3 );

   }else{
      return( 0 );

   }
}


/*====================================================================*\
   PolylineTransform() calculates the position of the input point in
   arc length - arc normal coordinates for the polyline specified by
   the input vertex list.  It returns the number of segments for which
   the point lies on a normal to the segment.

   The coordinates are calculated relative to all segments for which 
   input point lies on a normal to the segment.

   The line direction metric starts at the first point.  The line normal
   metric is positive on the right side of the line.

   Reginald H. Beardsley                            [e-mail]
\*====================================================================*/

#include <math.h>
#include <string.h>

int PolylineTransform( 
    float x0       /* X coordinate of input point */
   ,float y0       /* Y coordinate of input point */
   ,float* x       /* X coordinates of polyline vertices */
   ,float* y       /* Y coordinates of polyline vertices */
   ,float* s       /* arc length coordinate of polyline vertices */
   ,int n          /* number of vertices in polyline */
   ,float* il      /* arc length coordinate list for point */
   ,float* xl      /* crossline coordinate list for point */
   ){

   float a;   /* intercept of line */
   float b;   /* slope of line     */

   float x_;  /* X coordinate of intersection of normal & segment */
   float y_;  /* Y coordinate of intersection of normal & segment */

   float m_;  /* tangent of angle between m1 & m2                   */
   float m1;  /* vector from first point on segment to intersection */
   float m2;  /* vector from first point on segment to point        */

   float dx;  /* distance of point from line segement along normal  */

   int i;
   int nl = 0;

   memset( il ,0 ,sizeof(float)*n );
   memset( xl ,0 ,sizeof(float)*n );

   for( i=1; i<n; i++ ){

      if( x[i] != x[i-1] ){

         /*---------------------------*/
         /* non-vertical line segment */
         /*---------------------------*/

         b = (y[i] - y[i-1]) / (x[i] - x[i-1] );

         a = y[i-1] - x[i-1] * b;

         x_ = (x0 + (y0 - a) * b) / (1.0 + b*b);

         y_ = a + b * x_;

         m1 = (y_ - y[i-1]) / (x_ - x[i-1]);
         m2 = (y0 - y[i-1]) / (x0 - x[i-1]);
   
         m_ = (m1 - m2) / (1.0 + m1*m2);

         if( x[i] >= x_ && x_ >= x[i-1] 
          || x[i] <= x_ && x_ <= x[i-1] ){

            il[nl] = s[i-1] + sqrt( (y_-y[i-1])*(y_-y[i-1]) 
                                  + (x_-x[i-1])*(x_-x[i-1]) );

            dx = sqrt( (x0-x_)*(x0-x_) + (y0-y_)*(y0-y_) );

            if( m_ >= 0.0 ){
               xl[nl] = dx;

            }else{
               xl[nl] = dx * -1.0;

            }

            nl++;
         }

      }else{

         /*-----------------------*/
         /* vertical line segment */
         /*-----------------------*/

         if( y[i] >= y0 && y0 >= y[i-1] 
          || y[i] <= y0 && y0 <= y[i-1] ){
         
            x_ = x[i];
            y_ = y0;

            il[nl] = s[i-1] + fabs( y_ - y[i-1] );

            dx = fabs( x0 - x_);

            if( y0 >= y[i-1] && x0 > x_ ){
               xl[nl] = dx * -1.0;

            }else{
               xl[nl] = dx;

            }
            nl++;
         }
      }
   }

   return( nl );

}


void changeval(String type, Value *val, int f) {

	switch (*type) {
        case 's':
                err("can't change char header word");
        break;
        case 'h':
                val->h = f;
        break;
        case 'u':
                val->u = f;
        break;
        case 'l':
                val->l = f;
        break;
        case 'v':
                val->v = f;
        break;
        case 'i':
                val->i = f;
        break;
        case 'p':
                val->p = f;
        break;
        case 'f':
                val->f = f;
        break;
        case 'd':
                val->d = f;
        break;
        default:
                err("unknown type %s", type);
        break;
        }
}