Changeset 22456


Ignore:
Timestamp:
02/23/18 11:52:37 (7 years ago)
Author:
Mathieu Morlighem
Message:

CHG: cosmetics

File:
1 edited

Legend:

Unmodified
Added
Removed
  • issm/trunk-jpl/src/c/modules/ExpToLevelSetx/ExpToLevelSetxt.cpp

    r21931 r22456  
    99#endif
    1010
    11 /*Include files: {{{*/
     11/*Include files*/
    1212#include "./ExpToLevelSetx.h"
    1313double minimum_distance(double x1, double y1, double x2, double y2, double x0, double y0);
    14 void ContourToLevelSet(double* distance,double* contourx, double* contoury, int contournods, double* x, double* y, int i0, int i10);
    15 /*}}}*/
     14void   ContourToLevelSet(double* distance,double* contourx, double* contoury, int contournods, double* x, double* y, int i0, int i10);
    1615
    1716void* ExpToLevelSetxt(void* vpthread_handle){
     
    2928
    3029        /*recover parameters :*/
    31         Contours* contours  = gate->contours;
    32         int       nods      = gate->nods;
    33         double   *distance    = gate->distance;
    34         double   *x         = gate->x;
    35         double   *y         = gate->y;
     30        Contours* contours = gate->contours;
     31        int       nods     = gate->nods;
     32        double   *distance = gate->distance;
     33        double   *x        = gate->x;
     34        double   *y        = gate->y;
    3635
    3736        /*distribute indices across threads :*/
     
    3938
    4039        /*Loop through all contours: */
    41         for (i=0;i<contours->Size();i++){
     40        for(i=0;i<contours->Size();i++){
    4241                Contour<double>* contour=(Contour<double>*)contours->GetObjectByOffset(i);
    4342                ContourToLevelSet(distance,contour->x,contour->y,contour->nods,x,y,i0,i1);
     
    4847
    4948void ContourToLevelSet(double* dist,double* contourx, double* contoury, int contournods, double* x, double* y, int i0, int i1){/*{{{*/
    50         int i,j;
    5149        double x0,y0;
    5250        double x1,y1;
     
    5452        double mind;
    5553       
    56         for(i=i0;i<i1;i++){
     54        for(int i=i0;i<i1;i++){
     55
     56      /*Get current point*/
    5757                x0=x[i]; y0=y[i];
    5858
    5959                /*Figure out distance from (x0,y0) to contour: */
    6060                mind=1e+50;
    61                 for (j=0;j<contournods-1;j++){
    62                         x1=contourx[j]; y1=contoury[j];
     61                for(int j=0;j<contournods-1;j++){
     62         /*Get distance from current segment*/
     63                        x1=contourx[j];   y1=contoury[j];
    6364                        x2=contourx[j+1]; y2=contoury[j+1];
    6465                        mind=min(mind,minimum_distance(x1,y1,x2,y2,x0,y0));
     
    6667                dist[i]=min(dist[i],mind);
    6768        }
    68 
    6969}
    7070double ddistance(double x1,double y1,double x2,double y2){ return sqrt(pow(x2-x1,2)+pow(y2-y1,2)); }
    7171double ddot(double x1, double y1, double x2, double y2){ return x1*x2+y1*y2; }
    7272double minimum_distance(double x1, double y1, double x2, double y2, double x0, double y0){
     73        /* Return minimum distance between line segment [(x1,y1) (x2,y2)] and point (x0,y0) (v=(x1,y1), w=(x2,y2) and p=(x0,y0)*/
     74
     75   /*Get segment length square (avoid sqrt) |w-v|^2*/
     76        double l2 = pow(x2-x1,2)+pow(y2-y1,2);
     77
     78   /*segment is single point: v == w*/
     79        if(l2 == 0.) return ddistance(x0,y0,x1,y1);
     80
     81        /*Consider the line extending the segment, parameterized as v + t (w - v).
     82        We find projection of point p onto the line.
     83        It falls where t = [(p-v) . (w-v)] / |w-v|^2*/
     84        double t = ddot(x0-x1,y0-y1, x2-x1, y2-y1) / l2;
     85        if(t < 0.){
     86      /*Beyond the 'v' end of the segment*/
     87      return ddistance(x0,y0, x1, y1);
     88   }
     89        else if(t > 1.){
     90      /*Beyond the 'w' end of the segment*/
     91      return ddistance(x0,y0, x2,y2);
     92   }
    7393       
    74         // Return minimum distance between line segment [(x1,y1) (x2,y2)] and point (x0,y0) (v=(x1,y1), w=(x2,y2) and p=(x0,y0)
    75         double projectionx;
    76         double projectiony;
    77         double l2;
    78         double t;
    79 
    80         l2 = pow(x2-x1,2)+pow(y2-y1,2); // i.e. |w-v|^2 -  avoid a sqrt
    81 
    82         if (l2 == 0.0) return ddistance(x0,y0, x1,y1); // v == w case
    83         // Consider the line extending the segment, parameterized as v + t (w - v).
    84         //         // We find projection of point p onto the line.
    85         //           // It falls where t = [(p-v) . (w-v)] / |w-v|^2
    86         t = ddot(x0-x1,y0-y1, x2-x1, y2-y1) / l2;
    87         if (t < 0.0) return ddistance(x0,y0, x1, y1);       // Beyond the 'v' end of the segment
    88         else if (t > 1.0) return ddistance(x0,y0, x2,y2);  // Beyond the 'w' end of the segment
    89        
    90         projectionx= x1 + t* (x2-x1);  // Projection falls on the segment
    91         projectiony= y1 + t* (y2-y1);
     94   /*Projection falls on segment!*/
     95        double projectionx= x1 + t* (x2-x1);
     96        double projectiony= y1 + t* (y2-y1);
    9297        return ddistance(x0, y0, projectionx, projectiony);
    9398}
Note: See TracChangeset for help on using the changeset viewer.