Ice Sheet System Model  4.18
Code documentation
Functions
ExpSimplify.cpp File Reference
#include "./ExpSimplify.h"

Go to the source code of this file.

Functions

void ExpSimplifyUsage (void)
 
void simplify (Contour< double > *contour, bool *flags, int ind0, int ind1, double tolerance)
 
 WRAPPER (ExpSimplify_python)
 

Function Documentation

◆ ExpSimplifyUsage()

void ExpSimplifyUsage ( void  )

Definition at line 6 of file ExpSimplify.cpp.

6  {/*{{{*/
7  _printf_("ExpSimplify - Simplify Exp contour\n");
8  _printf_("\n");
9  _printf_(" Recursive Douglas-Peucker Polygon Simplification\n");
10  _printf_("\n");
11  _printf_(" Usage:\n");
12  _printf_(" ExpSimplify(expfile,tol);\n");
13  _printf_(" - expfile: name of the exp file\n");
14  _printf_(" - tol: tolerance (maximal euclidean distance allowed between the new line and a vertex)\n");
15  _printf_(" Additional options:\n");
16  _printf_(" - 'min': minimum number of vertices to save contours in exp file (default is 3)\n");
17  _printf_("\n");
18  _printf_(" Example:\n");
19  _printf_(" ExpSimplify('file.exp',100);\n");
20  _printf_(" ExpSimplify('file.exp',100,'min',4);\n");
21 }/*}}}*/

◆ simplify()

void simplify ( Contour< double > *  contour,
bool *  flags,
int  ind0,
int  ind1,
double  tolerance 
)

Definition at line 22 of file ExpSimplify.cpp.

22  {/*{{{*/
23 
24  bool closed = false;
25  double distance,beta,dx,dy;
26  double maxdistance;
27  int index = -1;
28  double *x = contour->x;
29  double *y = contour->y;
30 
31  /*Some checks*/
32  _assert_(ind0>=0 && ind0<contour->nods);
33  _assert_(ind1>=0 && ind1<contour->nods);
34  _assert_(ind1-ind0>=0);
35 
36  /*Check wether this portion is closed*/
37  if(x[ind0]==x[ind1] && y[ind0]==y[ind1]) closed=true;
38 
39  if(closed){
40 
41  /*calculate the shortest distance of all vertices between ind0 and ind1*/
42  for(int i=ind0;i<ind1-1;i++){
43  distance = sqrt((x[ind0]-x[i+1])*(x[ind0]-x[i+1]) + (y[ind0]-y[i+1])*(y[ind0]-y[i+1]));
44  if(i==ind0 || distance>maxdistance){
45  maxdistance=distance;
46  index = i + 1;
47  }
48  }
49  }
50  else{
51  /*calculate shortest distance of all points to the line from ind0 to ind1
52  * subtract starting point from other locations
53  *
54  * d = || (x-x0) - beta (xend - x0) ||
55  * <x-x0,xend-x0> = ||x-x0|| ||xend-x0|| cos(alpha)
56  * beta ||xend-x0|| = ||x-x0|| cos(alpha)
57  *
58  * So: beta = <x-x0,xend-x0>/<xend-x0,xend-x0> */
59 
60  for(int i=ind0+1;i<=ind1;i++){
61  beta = ((x[i]-x[ind0])*(x[ind1]-x[ind0]) + (y[i]-y[ind0])*(y[ind1]-y[ind0]))/((x[ind1]-x[ind0])*(x[ind1]-x[ind0])+(y[ind1]-y[ind0])*(y[ind1]-y[ind0]));
62  dx = x[i]-beta*x[ind1]+(beta-1.)*x[ind0];
63  dy = y[i]-beta*y[ind1]+(beta-1.)*y[ind0];
64  distance = sqrt(dx*dx + dy*dy);
65  if(i==ind0+1 || distance>maxdistance){
66  maxdistance = distance;
67  index = i;
68  }
69  }
70 
71  }
72 
73  /*if the maximum distance is smaller than the tolerance remove vertices between ind0 and ind1*/
74  if(maxdistance<tolerance){
75  if(ind0!=ind1-1){
76  for(int i=ind0+1;i<ind1;i++) flags[i]=false;
77  }
78  }
79  else{
80  /*if not, call simplifyrec for the segments between ind0 and index
81  * (index and ind1)*/
82  _assert_(index!=-1);
83  _assert_(index!=ind1);
84  _assert_(index!=ind0);
85  simplify(contour,flags,ind0 ,index,tolerance);
86  simplify(contour,flags,index,ind1, tolerance);
87  }
88 
89 }/*}}}*/

◆ WRAPPER()

WRAPPER ( ExpSimplify_python  )

Definition at line 90 of file ExpSimplify.cpp.

90  {
91 
92  int i,verbose=1;
93 
94  /*input: */
95  char* expfile = NULL;
96  double tolerance;
97  Options *options = NULL;
98  double minimumvertices_double;
99  int minimumvertices;
100 
101  /*output*/
102  Contours* oldcontours = NULL;
103  Contours* newcontours = NULL;
104 
105  /*Boot module: */
106  MODULEBOOT();
107 
108  /*checks on arguments: */
109  /*checks on arguments on the matlab side: */
110  if (nrhs<NRHS || nlhs>NLHS){
111  ExpSimplifyUsage(); _error_("ExpSimplify usage error");
112  }
113 
114  /*Input datasets: */
115  FetchData(&expfile, EXPFILE);
116  FetchData(&tolerance,TOLERANCE);
117  FetchData(&options,NRHS,nrhs,prhs);
118 
119  /*some checks*/
120  if(tolerance<0) _error_("tolerance must be a positivve scalar");
121 
122  /* Run core computations: */
123  Contour<double>* contour = NULL;
124  Contour<double>* newcontour = NULL;
125  int nods,newnods;
126  bool* flags = NULL;
127  double* x = NULL;
128  double* y = NULL;
129  double distance;
130 
131  /*Process options*/
132  options->Get(&minimumvertices_double,"min",3.);
133  if(minimumvertices_double<1.) _error_("'min' (minimum number of verties) should be a positive integer");
134  minimumvertices = int(minimumvertices_double);
135 
136  /*Read old contours and allocate new contours*/
137  oldcontours=ExpRead<double>(expfile);
138  newcontours=new Contours();
139  for(int counter=0;counter<oldcontours->Size();counter++){
140 
141  /*Get single contour*/
142  contour = (Contour<double>*)oldcontours->GetObjectByOffset(counter);
143  nods = contour->nods;
144  x = contour->x;
145  y = contour->y;
146  _printf_(" Initial number of vertices in contour #"<<counter+1<<": "<<nods << "\n");
147 
148  /*Allocate flags (1=keep, 0=remove)*/
149  if(nods>0) flags = xNew<bool>(nods);
150 
151  if(nods==0){
152  /*Don't do anything*/
153  }
154  else if(nods==1){
155  flags[0] = true;
156  }
157  else if(nods==2){
158  /*check if the distance between both is less than the tolerance
159  * If so, return the center*/
160  distance = sqrt((x[1]-x[0])*(x[1]-x[0]) + (y[1]-y[0])*(y[1]-y[0]));
161  if(distance<=tolerance){
162  x[0]=(x[0]+x[1])/2.;
163  y[0]=(y[0]+y[1])/2.;
164  flags[0] = true;
165  flags[1] = false;
166  }
167  else{
168  flags[0] = true;
169  flags[1] = true;
170  }
171  }
172  else{
173  /*Start recursive call to simplify*/
174  for(int i=0;i<nods;i++) flags[i]=true;
175  simplify(contour,flags,0,nods-1,tolerance);
176  }
177 
178  /*Add new contour to newcontours*/
179  newnods = 0;
180  for(int i=0;i<nods;i++) if(flags[i]) newnods++;
181 
182  /*Do we save new profile?*/
183  if(newnods>=minimumvertices){
184  _printf_(" Final number of vertices in contour #"<<counter+1<<": "<<newnods << "\n");
185  newcontour = new Contour<double>();
186  newcontour->nods = newnods;
187  newcontour->x = xNew<double>(newnods);
188  newcontour->y = xNew<double>(newnods);
189  newnods = 0;
190  for(int i=0;i<nods;i++){
191  if(flags[i]){
192  newcontour->x[newnods] = contour->x[i];
193  newcontour->y[newnods] = contour->y[i];
194  newnods++;
195  }
196  }
197  _assert_(newnods==newcontour->nods);
198 
199  /*Add to main dataset*/
200  newcontours->AddObject(newcontour);
201  }
202  else{
203  _printf_(" Final number of vertices in contour #"<<counter+1<<": "<<newnods<<" (not saved)\n");
204  }
205 
206  /*cleanup*/
207  xDelete<bool>(flags);
208  }
209  _printf_(" Initial number of contours: "<<oldcontours->Size() << "\n");
210  _printf_(" Final number of contours: "<<newcontours->Size() << "\n");
211 
212  /*Write data: */
213  ExpWrite(newcontours,expfile);
214 
215  /*Clean-up*/
216  xDelete<char>(expfile);
217  delete options;
218  delete oldcontours;
219  delete newcontours;
220 
221  /*end module: */
222  MODULEEND();
223 }
DataSet::Size
int Size()
Definition: DataSet.cpp:399
Options
Definition: Options.h:9
_assert_
#define _assert_(ignore)
Definition: exceptions.h:37
DataSet::AddObject
int AddObject(Object *object)
Definition: DataSet.cpp:252
_printf_
#define _printf_(StreamArgs)
Definition: Print.h:22
NRHS
#define NRHS
Definition: BamgConvertMesh.h:52
Contours
Declaration of Contours class.
Definition: Contours.h:10
Contour
Definition: Contour.h:15
ExpWrite
int ExpWrite(Contours *contours, char *domainname)
Definition: Contours.cpp:32
FetchData
void FetchData(char **pstring, char *stringin)
Definition: FetchJavascriptData.cpp:16
simplify
void simplify(Contour< double > *contour, bool *flags, int ind0, int ind1, double tolerance)
Definition: ExpSimplify.cpp:22
Contour::nods
int nods
Definition: Contour.h:20
Options::Get
void Get(OptionType *pvalue, const char *name)
Definition: Options.h:21
NLHS
#define NLHS
Definition: BamgConvertMesh.h:50
Contour::y
doubletype * y
Definition: Contour.h:22
_error_
#define _error_(StreamArgs)
Definition: exceptions.h:49
DataSet::GetObjectByOffset
Object * GetObjectByOffset(int offset)
Definition: DataSet.cpp:334
Contour::x
doubletype * x
Definition: Contour.h:21
ExpSimplifyUsage
void ExpSimplifyUsage(void)
Definition: ExpSimplify.cpp:6