Statistics
| Revision:

root / org.gvsig.toolbox / trunk / org.gvsig.toolbox / org.gvsig.toolbox.algorithm / src / main / java / es / unex / sextante / tin / smoothTinBezier / Bezier2.java @ 59

History | View | Annotate | Download (36.1 KB)

1
/**
2
 *    @author              Josef Bezdek, ZCU Plzen
3
 *          @version             1.0
4
 *    @since                 JDK1.5
5
 */
6

    
7
package es.unex.sextante.tin.smoothTinBezier;
8

    
9
import java.util.Iterator;
10
import java.util.LinkedList;
11

    
12
import com.vividsolutions.jts.geom.Coordinate;
13
import com.vividsolutions.jts.geom.Envelope;
14
import com.vividsolutions.jts.geom.GeometryFactory;
15
import com.vividsolutions.jts.geom.LinearRing;
16
import com.vividsolutions.jts.geom.impl.CoordinateArraySequence;
17

    
18
class Bezier2 {
19
   Coordinate normalN1 = null;
20
   Coordinate normalN2 = null;
21
   Coordinate normalN3 = null;
22
   Coordinate b300;
23
   Coordinate b030;
24
   Coordinate b003;
25
   Coordinate b012     = null;
26
   Coordinate b021     = null;
27
   Coordinate b102     = null;
28
   Coordinate b120     = null;
29
   Coordinate b210     = null;
30
   Coordinate b201     = null;
31
   Coordinate A111     = null;
32
   Coordinate B111     = null;
33
   Coordinate C111     = null;
34

    
35
   Coordinate G        = null;
36

    
37
   Coordinate A201, A102, A012, A021, B201, B102;
38

    
39
   Coordinate n200     = null;
40
   Coordinate n020     = null;
41
   Coordinate n002     = null;
42
   Coordinate n110     = null;
43
   Coordinate n011     = null;
44
   Coordinate n101     = null;
45

    
46

    
47
   /******************************************************************************************************************************
48
    * Constructor
49
    * 
50
    * @param T
51
    *                Trinagle fom original TIN
52
    * @param listA -
53
    *                normal vectors of planes in point A of triangle T
54
    * @param listB -
55
    *                normal vectors of planes in point B of triangle T
56
    * @param listC -
57
    *                normal vectors of planes in point C of triangle T
58
    */
59
   Bezier2(final TriangleDT T,
60
           final LinkedList listA,
61
           final LinkedList listB,
62
           final LinkedList listC,
63
           final int typeOfBreakLine) {
64
      b300 = T.A;
65
      b030 = T.B;
66
      b003 = T.C;
67
      setNormalVector(listA, listB, listC);
68
      setControlPoints(typeOfBreakLine);
69

    
70
   }
71

    
72

    
73
   /******************************************************************************************************************************
74
    * Constructor
75
    * 
76
    * @param coords -
77
    *                array of three vertexes which generate triangle
78
    */
79
   Bezier2(final Coordinate[] coords) {
80
      b300 = coords[0];
81
      b030 = coords[1];
82
      b003 = coords[2];
83
   }
84

    
85

    
86
   /******************************************************************************************************************************
87
    * Constructor
88
    * 
89
    * @param coords -
90
    *                array of three vertexes which generate triangle
91
    * @param listA -
92
    *                list of normal vektors around vertex A
93
    * @param listB -
94
    *                list of normal vektors around vertex B
95
    * @param listC -
96
    *                list of normal vektors around vertex C
97
    * @param typeOfBreakLine -
98
    *                type of break line which triangle contains
99
    */
100
   Bezier2(final Coordinate[] coords,
101
           final LinkedList listA,
102
           final LinkedList listB,
103
           final LinkedList listC,
104
           final int typeOfBreakLine) {
105
      b300 = coords[0];
106
      b030 = coords[1];
107
      b003 = coords[2];
108
      setNormalVector(listA, listB, listC);
109
      setControlPoints(typeOfBreakLine);
110
   }
111

    
112

    
113
   /******************************************************************************************************************************
114
    * The Protected method for setting one normals for each vertex of T
115
    * 
116
    * @param listA -
117
    *                normal vectors of planes in point A of triangle T
118
    * @param listB -
119
    *                normal vectors of planes in point B of triangle T
120
    * @param listC -
121
    *                normal vectors of planes in point C of triangle T
122
    */
123
   protected void setNormalVector(final LinkedList listA,
124
                                  final LinkedList listB,
125
                                  final LinkedList listC) {
126
      normalN1 = countVector(listA, b300);
127
      normalN2 = countVector(listB, b030);
128
      normalN3 = countVector(listC, b003);
129
   }
130

    
131

    
132
   /******************************************************************************************************************************
133
    * The method which sets new vecter between points A and B
134
    * 
135
    * @param A -
136
    *                start point
137
    * @param B -
138
    *                stop point
139
    * @return vector AB
140
    */
141
   protected static Coordinate setVector(final Coordinate A,
142
                                         final Coordinate B) {
143
      return new Coordinate(B.x - A.x, B.y - A.y, B.z - A.z);
144

    
145
   }
146

    
147

    
148
   /******************************************************************************************************************************
149
    * The method for setting normal vectors of two vectors
150
    * 
151
    * @param A -
152
    *                vector A
153
    * @param B -
154
    *                vector B
155
    * @return normal vector
156
    */
157
   protected static Coordinate setNormalVector(final Coordinate A,
158
                                               final Coordinate B) {
159
      final Coordinate normal = new Coordinate(A.y * B.z - A.z * B.y, A.z * B.x - A.x * B.z, (A.x * B.y - A.y * B.x));
160
      final double sum = Math.sqrt(Math.pow(normal.x, 2) + Math.pow(normal.y, 2) + Math.pow(normal.z, 2));
161
      //double sum = 1;
162
      if (normal.z > 0) {
163
         return new Coordinate((normal.x / sum), (normal.y / sum), (normal.z / sum));
164
      }
165
      else {
166
         return new Coordinate((-1) * (normal.x / sum), (-1) * (normal.y / sum), (-1) * (normal.z / sum));
167
      }
168
   }
169

    
170

    
171
   /******************************************************************************************************************************
172
    * The private method counts one normal vector from normal vectors of every plane in vertex of triangle
173
    * 
174
    * @param list -
175
    *                normal vectors of planes in point of triangle T
176
    * @param P /
177
    *                vertex of T
178
    * @return normal vector
179
    */
180
   private static Coordinate countVector(final LinkedList list,
181
                                         final Coordinate P) {
182
      final Iterator iter = list.iterator();
183
      //double koeficient = 1D;
184
      double sumX = 0;
185
      double sumY = 0;
186
      double sumZ = 0;
187
      while (iter.hasNext()) {
188
         final Coordinate X = (Coordinate) iter.next();
189
         sumX += X.x;
190
         sumY += X.y;
191
         sumZ += X.z;
192
      }
193
      final double sum = Math.sqrt(Math.pow(sumX, 2) + Math.pow(sumY, 2) + Math.pow(sumZ, 2));
194
      return new Coordinate((sumX / sum), (sumY / sum), (sumZ / sum));
195
   }
196

    
197

    
198
   /******************************************************************************************************************************
199
    * The protected method counts Scalar product of two vectors v1,v2
200
    * 
201
    * @param v1 -
202
    *                vector
203
    * @param v2 -
204
    *                vector
205
    * @return - scalar product
206
    */
207
   protected static double countScalarProduct(final Coordinate v1,
208
                                              final Coordinate v2) {
209
      final double scalar = v1.x * v2.x + v1.y * v2.y + v1.z * v2.z;
210
      return scalar;
211
   }
212

    
213

    
214
   /******************************************************************************************************************************
215
    * The protected function counts cross vector product of two vectors A,B
216
    * 
217
    * @param A -
218
    *                vector
219
    * @param B -
220
    *                vector
221
    * @return - cross product
222
    */
223
   protected static Coordinate countCrossProduct(final Coordinate A,
224
                                                 final Coordinate B) {
225
      final Coordinate normal = new Coordinate(A.y * B.z - A.z * B.y, A.z * B.x - A.x * B.z, (A.x * B.y - A.y * B.x));
226
      final double sum = Math.sqrt(Math.pow(normal.x, 2) + Math.pow(normal.y, 2) + Math.pow(normal.z, 2));
227
      if (normal.z > 0) {
228
         return new Coordinate((normal.x / sum), (normal.y / sum), (normal.z / sum));
229
      }
230
      else {
231
         return new Coordinate((-1) * (normal.x / sum), (-1) * (normal.y / sum), (-1) * (normal.z / sum));
232
      }
233

    
234
   }
235

    
236

    
237
   /******************************************************************************************************************************
238
    * The protected method counts difference of two vectors v1, v2
239
    * 
240
    * @param v1 -
241
    *                vector
242
    * @param v2 -
243
    *                vector
244
    * @return difference vector
245
    */
246
   protected static Coordinate countDifferenceProduct(final Coordinate v1,
247
                                                      final Coordinate v2) {
248
      return new Coordinate(v1.x - v2.x, v1.y - v2.y, v1.z - v2.z);
249
   }
250

    
251

    
252
   /******************************************************************************************************************************
253
    * The protected function counts sum of two vectors v1, v2
254
    * 
255
    * @param v1 -
256
    *                vector
257
    * @param v2 -
258
    *                vector
259
    * @return sum vector
260
    */
261
   protected static Coordinate countSumProduct(final Coordinate v1,
262
                                               final Coordinate v2) {
263
      return new Coordinate(v1.x + v2.x, v1.y + v2.y, v1.z + v2.z);
264
   }
265

    
266

    
267
   /******************************************************************************************************************************
268
    * The protected function which normalize vector v
269
    * 
270
    * @param v -
271
    *                vector
272
    * @return normalized vector
273
    */
274
   protected static Coordinate normalizeVect(final Coordinate v) {
275
      final double sum = Math.sqrt(Math.pow(v.x, 2) + Math.pow(v.y, 2) + Math.pow(v.z, 2));
276
      return new Coordinate((v.x / sum), (v.y / sum), (v.z / sum));
277
   }
278

    
279

    
280
   /******************************************************************************************************************************
281
    * The private function for help calculation
282
    */
283
   private double helpCount(final Coordinate Pi,
284
                            final Coordinate Pj,
285
                            final Coordinate Ni,
286
                            final Coordinate Nj) {
287
      return 2 * (countScalarProduct(countDifferenceProduct(Pj, Pi), countSumProduct(Ni, Nj)) / countScalarProduct(
288
               countDifferenceProduct(Pj, Pi), countDifferenceProduct(Pj, Pi)));
289
   }
290

    
291

    
292
   /******************************************************************************************************************************
293
    * The protected method sets Quadratic normals of Bezier triangle
294
    */
295
   protected void setQuadraticNormals() {
296
      n200 = normalN1;
297
      n020 = normalN2;
298
      n002 = normalN3;
299

    
300

    
301
      double help = helpCount(b300, b030, normalN1, normalN2);
302
      Coordinate dif = countDifferenceProduct(b030, b300);
303
      dif.x = dif.x * help;
304
      dif.y = dif.y * help;
305
      dif.z = dif.z * help;
306
      n110 = normalizeVect(countDifferenceProduct(countSumProduct(normalN1, normalN2), dif));
307
      help = helpCount(b030, b003, normalN2, normalN3);
308
      dif = countDifferenceProduct(b003, b030);
309
      dif.x = dif.x * help;
310
      dif.y = dif.y * help;
311
      dif.z = dif.z * help;
312
      n011 = normalizeVect(countDifferenceProduct(countSumProduct(normalN2, normalN3), dif));
313
      help = helpCount(b003, b300, normalN3, normalN1);
314
      dif = countDifferenceProduct(b300, b003);
315
      dif.x = dif.x * help;
316
      dif.y = dif.y * help;
317
      dif.z = dif.z * help;
318
      n101 = normalizeVect(countDifferenceProduct(countSumProduct(normalN3, normalN1), dif));
319
   }
320

    
321

    
322
   /******************************************************************************************************************************
323
    * The protected method gets normal into Bezier triangle with barycentric coordinate
324
    * 
325
    * @param u -
326
    *                barycentric coordinate u
327
    * @param v -
328
    *                barycentric coordinate u
329
    * @return normal of surface
330
    */
331
   protected Coordinate getNormal(final double u,
332
                                  final double v) {
333
      final double w = 1 - (u + v);
334
      final double x = n200.x * Math.pow(w, 2) + n020.x * Math.pow(u, 2) + n002.x * Math.pow(v, 2) + n110.x * w * u + n011.x * u
335
                       * v + n101.x * w * v;
336

    
337
      final double y = n200.y * Math.pow(w, 2) + n020.y * Math.pow(u, 2) + n002.y * Math.pow(v, 2) + n110.y * w * u + n011.y * u
338
                       * v + n101.y * w * v;
339

    
340
      final double z = n200.z * Math.pow(w, 2) + n020.z * Math.pow(u, 2) + n002.z * Math.pow(v, 2) + n110.z * w * u + n011.z * u
341
                       * v + n101.z * w * v;
342
      return new Coordinate(x, y, z);
343
   }
344

    
345

    
346
   /******************************************************************************************************************************
347
    * The protected method counts project point to plane
348
    * 
349
    * @param pointOfPlane -
350
    *                point which is contained into plane
351
    * @param normalOfPlane -
352
    *                normal vector of plane
353
    * @param pointOfLine -
354
    *                point which will project to plane
355
    * @param normalOfLine -
356
    *                normal wich get direction of projection to plane
357
    * @return coordinates of projected point in plane
358
    */
359
   protected Coordinate countProjectOnToPlane(final Coordinate pointOfPlane,
360
                                              final Coordinate normalOfPlane,
361
                                              final Coordinate pointOfLine,
362
                                              final Coordinate normalOfLine) {
363
      final double d = -normalOfPlane.x * pointOfPlane.x - normalOfPlane.y * pointOfPlane.y - normalOfPlane.z * pointOfPlane.z;
364
      final double param = (-pointOfLine.x * normalOfPlane.x - pointOfLine.y * normalOfPlane.y - pointOfLine.z * normalOfPlane.z - d)
365
                           / (normalOfLine.x * normalOfPlane.x + normalOfLine.y * normalOfPlane.y + normalOfLine.z
366
                                                                                                    * normalOfPlane.z);
367
      return new Coordinate(pointOfLine.x + param * normalOfLine.x, pointOfLine.y + param * normalOfLine.y, pointOfLine.z + param
368
                                                                                                            * normalOfLine.z);
369

    
370
   }
371

    
372

    
373
   /******************************************************************************************************************************
374
    * The method for setting control points of bezier triangle
375
    * 
376
    * @param typeOfBreakLine -
377
    *                type of break which triangle contains
378
    */
379
   protected void setControlPoints(final int typeOfBreakLine) {
380
      setQuadraticNormals();
381

    
382
      switch (typeOfBreakLine) {
383
         case (-1): {
384
            b210 = new Coordinate((2 * b300.x + b030.x - countScalarProduct(countDifferenceProduct(b030, b300), normalN1)
385
                                                         * normalN1.x) / 3, (2 * b300.y + b030.y - countScalarProduct(
386
                     countDifferenceProduct(b030, b300), normalN1)
387
                                                                                                   * normalN1.y) / 3,
388
                     (2 * b300.z + b030.z - countScalarProduct(countDifferenceProduct(b030, b300), normalN1) * normalN1.z) / 3);
389

    
390
            b120 = new Coordinate((2 * b030.x + b300.x - countScalarProduct(countDifferenceProduct(b300, b030), normalN2)
391
                                                         * normalN2.x) / 3, (2 * b030.y + b300.y - countScalarProduct(
392
                     countDifferenceProduct(b300, b030), normalN2)
393
                                                                                                   * normalN2.y) / 3,
394
                     (2 * b030.z + b300.z - countScalarProduct(countDifferenceProduct(b300, b030), normalN2) * normalN2.z) / 3);
395

    
396
            b021 = new Coordinate((2 * b030.x + b003.x - countScalarProduct(countDifferenceProduct(b003, b030), normalN2)
397
                                                         * normalN2.x) / 3, (2 * b030.y + b003.y - countScalarProduct(
398
                     countDifferenceProduct(b003, b030), normalN2)
399
                                                                                                   * normalN2.y) / 3,
400
                     (2 * b030.z + b003.z - countScalarProduct(countDifferenceProduct(b003, b030), normalN2) * normalN2.z) / 3);
401

    
402
            b012 = new Coordinate((2 * b003.x + b030.x - countScalarProduct(countDifferenceProduct(b030, b003), normalN3)
403
                                                         * normalN3.x) / 3, (2 * b003.y + b030.y - countScalarProduct(
404
                     countDifferenceProduct(b030, b003), normalN3)
405
                                                                                                   * normalN3.y) / 3,
406
                     (2 * b003.z + b030.z - countScalarProduct(countDifferenceProduct(b030, b003), normalN3) * normalN3.z) / 3);
407

    
408
            b201 = new Coordinate((2 * b300.x + b003.x - countScalarProduct(countDifferenceProduct(b003, b300), normalN1)
409
                                                         * normalN1.x) / 3, (2 * b300.y + b003.y - countScalarProduct(
410
                     countDifferenceProduct(b003, b300), normalN1)
411
                                                                                                   * normalN1.y) / 3,
412
                     (2 * b300.z + b003.z - countScalarProduct(countDifferenceProduct(b003, b300), normalN1) * normalN1.z) / 3);
413

    
414
            b102 = new Coordinate((2 * b003.x + b300.x - countScalarProduct(countDifferenceProduct(b300, b003), normalN3)
415
                                                         * normalN3.x) / 3, (2 * b003.y + b300.y - countScalarProduct(
416
                     countDifferenceProduct(b300, b003), normalN3)
417
                                                                                                   * normalN3.y) / 3,
418
                     (2 * b003.z + b300.z - countScalarProduct(countDifferenceProduct(b300, b003), normalN3) * normalN3.z) / 3);
419
            break;
420
         }
421
         case (0): {
422
            b210 = new Coordinate((2 * b300.x + b030.x) / 3, (2 * b300.y + b030.y) / 3, (2 * b300.z + b030.z) / 3);
423

    
424
            b120 = new Coordinate((2 * b030.x + b300.x) / 3, (2 * b030.y + b300.y) / 3, (2 * b030.z + b300.z) / 3);
425

    
426
            b021 = new Coordinate((2 * b030.x + b003.x - countScalarProduct(countDifferenceProduct(b003, b030), normalN2)
427
                                                         * normalN2.x) / 3, (2 * b030.y + b003.y - countScalarProduct(
428
                     countDifferenceProduct(b003, b030), normalN2)
429
                                                                                                   * normalN2.y) / 3,
430
                     (2 * b030.z + b003.z - countScalarProduct(countDifferenceProduct(b003, b030), normalN2) * normalN2.z) / 3);
431

    
432
            b012 = new Coordinate((2 * b003.x + b030.x - countScalarProduct(countDifferenceProduct(b030, b003), normalN3)
433
                                                         * normalN3.x) / 3, (2 * b003.y + b030.y - countScalarProduct(
434
                     countDifferenceProduct(b030, b003), normalN3)
435
                                                                                                   * normalN3.y) / 3,
436
                     (2 * b003.z + b030.z - countScalarProduct(countDifferenceProduct(b030, b003), normalN3) * normalN3.z) / 3);
437

    
438
            b201 = new Coordinate((2 * b300.x + b003.x - countScalarProduct(countDifferenceProduct(b003, b300), normalN1)
439
                                                         * normalN1.x) / 3, (2 * b300.y + b003.y - countScalarProduct(
440
                     countDifferenceProduct(b003, b300), normalN1)
441
                                                                                                   * normalN1.y) / 3,
442
                     (2 * b300.z + b003.z - countScalarProduct(countDifferenceProduct(b003, b300), normalN1) * normalN1.z) / 3);
443

    
444
            b102 = new Coordinate((2 * b003.x + b300.x - countScalarProduct(countDifferenceProduct(b300, b003), normalN3)
445
                                                         * normalN3.x) / 3, (2 * b003.y + b300.y - countScalarProduct(
446
                     countDifferenceProduct(b300, b003), normalN3)
447
                                                                                                   * normalN3.y) / 3,
448
                     (2 * b003.z + b300.z - countScalarProduct(countDifferenceProduct(b300, b003), normalN3) * normalN3.z) / 3);
449
            break;
450
         }
451
         case (1): {
452
            b210 = new Coordinate((2 * b300.x + b030.x - countScalarProduct(countDifferenceProduct(b030, b300), normalN1)
453
                                                         * normalN1.x) / 3, (2 * b300.y + b030.y - countScalarProduct(
454
                     countDifferenceProduct(b030, b300), normalN1)
455
                                                                                                   * normalN1.y) / 3,
456
                     (2 * b300.z + b030.z - countScalarProduct(countDifferenceProduct(b030, b300), normalN1) * normalN1.z) / 3);
457

    
458
            b120 = new Coordinate((2 * b030.x + b300.x - countScalarProduct(countDifferenceProduct(b300, b030), normalN2)
459
                                                         * normalN2.x) / 3, (2 * b030.y + b300.y - countScalarProduct(
460
                     countDifferenceProduct(b300, b030), normalN2)
461
                                                                                                   * normalN2.y) / 3,
462
                     (2 * b030.z + b300.z - countScalarProduct(countDifferenceProduct(b300, b030), normalN2) * normalN2.z) / 3);
463

    
464
            b021 = new Coordinate((2 * b030.x + b003.x) / 3, (2 * b030.y + b003.y) / 3, (2 * b030.z + b003.z) / 3);
465

    
466
            b012 = new Coordinate((2 * b003.x + b030.x) / 3, (2 * b003.y + b030.y) / 3, (2 * b003.z + b030.z) / 3);
467

    
468
            b201 = new Coordinate((2 * b300.x + b003.x - countScalarProduct(countDifferenceProduct(b003, b300), normalN1)
469
                                                         * normalN1.x) / 3, (2 * b300.y + b003.y - countScalarProduct(
470
                     countDifferenceProduct(b003, b300), normalN1)
471
                                                                                                   * normalN1.y) / 3,
472
                     (2 * b300.z + b003.z - countScalarProduct(countDifferenceProduct(b003, b300), normalN1) * normalN1.z) / 3);
473

    
474
            b102 = new Coordinate((2 * b003.x + b300.x - countScalarProduct(countDifferenceProduct(b300, b003), normalN3)
475
                                                         * normalN3.x) / 3, (2 * b003.y + b300.y - countScalarProduct(
476
                     countDifferenceProduct(b300, b003), normalN3)
477
                                                                                                   * normalN3.y) / 3,
478
                     (2 * b003.z + b300.z - countScalarProduct(countDifferenceProduct(b300, b003), normalN3) * normalN3.z) / 3);
479
            break;
480
         }
481
         case (2): {
482
            b210 = new Coordinate((2 * b300.x + b030.x - countScalarProduct(countDifferenceProduct(b030, b300), normalN1)
483
                                                         * normalN1.x) / 3, (2 * b300.y + b030.y - countScalarProduct(
484
                     countDifferenceProduct(b030, b300), normalN1)
485
                                                                                                   * normalN1.y) / 3,
486
                     (2 * b300.z + b030.z - countScalarProduct(countDifferenceProduct(b030, b300), normalN1) * normalN1.z) / 3);
487

    
488
            b120 = new Coordinate((2 * b030.x + b300.x - countScalarProduct(countDifferenceProduct(b300, b030), normalN2)
489
                                                         * normalN2.x) / 3, (2 * b030.y + b300.y - countScalarProduct(
490
                     countDifferenceProduct(b300, b030), normalN2)
491
                                                                                                   * normalN2.y) / 3,
492
                     (2 * b030.z + b300.z - countScalarProduct(countDifferenceProduct(b300, b030), normalN2) * normalN2.z) / 3);
493

    
494
            b021 = new Coordinate((2 * b030.x + b003.x - countScalarProduct(countDifferenceProduct(b003, b030), normalN2)
495
                                                         * normalN2.x) / 3, (2 * b030.y + b003.y - countScalarProduct(
496
                     countDifferenceProduct(b003, b030), normalN2)
497
                                                                                                   * normalN2.y) / 3,
498
                     (2 * b030.z + b003.z - countScalarProduct(countDifferenceProduct(b003, b030), normalN2) * normalN2.z) / 3);
499

    
500
            b012 = new Coordinate((2 * b003.x + b030.x - countScalarProduct(countDifferenceProduct(b030, b003), normalN3)
501
                                                         * normalN3.x) / 3, (2 * b003.y + b030.y - countScalarProduct(
502
                     countDifferenceProduct(b030, b003), normalN3)
503
                                                                                                   * normalN3.y) / 3,
504
                     (2 * b003.z + b030.z - countScalarProduct(countDifferenceProduct(b030, b003), normalN3) * normalN3.z) / 3);
505

    
506
            b201 = new Coordinate((2 * b300.x + b003.x) / 3, (2 * b300.y + b003.y) / 3, (2 * b300.z + b003.z) / 3);
507

    
508
            b102 = new Coordinate((2 * b003.x + b300.x) / 3, (2 * b003.y + b300.y) / 3, (2 * b003.z + b300.z) / 3);
509
            break;
510
         }
511
         case (3): {
512
            b210 = new Coordinate((2 * b300.x + b030.x) / 3, (2 * b300.y + b030.y) / 3, (2 * b300.z + b030.z) / 3);
513

    
514
            b120 = new Coordinate((2 * b030.x + b300.x) / 3, (2 * b030.y + b300.y) / 3, (2 * b030.z + b300.z) / 3);
515

    
516
            b021 = new Coordinate((2 * b030.x + b003.x - countScalarProduct(countDifferenceProduct(b003, b030), normalN2)
517
                                                         * normalN2.x) / 3, (2 * b030.y + b003.y - countScalarProduct(
518
                     countDifferenceProduct(b003, b030), normalN2)
519
                                                                                                   * normalN2.y) / 3,
520
                     (2 * b030.z + b003.z - countScalarProduct(countDifferenceProduct(b003, b030), normalN2) * normalN2.z) / 3);
521

    
522
            b012 = new Coordinate((2 * b003.x + b030.x - countScalarProduct(countDifferenceProduct(b030, b003), normalN3)
523
                                                         * normalN3.x) / 3, (2 * b003.y + b030.y - countScalarProduct(
524
                     countDifferenceProduct(b030, b003), normalN3)
525
                                                                                                   * normalN3.y) / 3,
526
                     (2 * b003.z + b030.z - countScalarProduct(countDifferenceProduct(b030, b003), normalN3) * normalN3.z) / 3);
527

    
528
            b201 = new Coordinate((2 * b300.x + b003.x) / 3, (2 * b300.y + b003.y) / 3, (2 * b300.z + b003.z) / 3);
529

    
530
            b102 = new Coordinate((2 * b003.x + b300.x) / 3, (2 * b003.y + b300.y) / 3, (2 * b003.z + b300.z) / 3);
531
            break;
532
         }
533

    
534
         case (5): {
535
            b210 = new Coordinate((2 * b300.x + b030.x - countScalarProduct(countDifferenceProduct(b030, b300), normalN1)
536
                                                         * normalN1.x) / 3, (2 * b300.y + b030.y - countScalarProduct(
537
                     countDifferenceProduct(b030, b300), normalN1)
538
                                                                                                   * normalN1.y) / 3,
539
                     (2 * b300.z + b030.z - countScalarProduct(countDifferenceProduct(b030, b300), normalN1) * normalN1.z) / 3);
540

    
541
            b120 = new Coordinate((2 * b030.x + b300.x - countScalarProduct(countDifferenceProduct(b300, b030), normalN2)
542
                                                         * normalN2.x) / 3, (2 * b030.y + b300.y - countScalarProduct(
543
                     countDifferenceProduct(b300, b030), normalN2)
544
                                                                                                   * normalN2.y) / 3,
545
                     (2 * b030.z + b300.z - countScalarProduct(countDifferenceProduct(b300, b030), normalN2) * normalN2.z) / 3);
546

    
547
            b021 = new Coordinate((2 * b030.x + b003.x) / 3, (2 * b030.y + b003.y) / 3, (2 * b030.z + b003.z) / 3);
548

    
549
            b012 = new Coordinate((2 * b003.x + b030.x) / 3, (2 * b003.y + b030.y) / 3, (2 * b003.z + b030.z) / 3);
550

    
551
            b201 = new Coordinate((2 * b300.x + b003.x) / 3, (2 * b300.y + b003.y) / 3, (2 * b300.z + b003.z) / 3);
552

    
553
            b102 = new Coordinate((2 * b003.x + b300.x) / 3, (2 * b003.y + b300.y) / 3, (2 * b003.z + b300.z) / 3);
554
            break;
555
         }
556
         case (4): {
557
            b210 = new Coordinate((2 * b300.x + b030.x) / 3, (2 * b300.y + b030.y) / 3, (2 * b300.z + b030.z) / 3);
558

    
559
            b120 = new Coordinate((2 * b030.x + b300.x) / 3, (2 * b030.y + b300.y) / 3, (2 * b030.z + b300.z) / 3);
560

    
561
            b021 = new Coordinate((2 * b030.x + b003.x) / 3, (2 * b030.y + b003.y) / 3, (2 * b030.z + b003.z) / 3);
562

    
563
            b012 = new Coordinate((2 * b003.x + b030.x) / 3, (2 * b003.y + b030.y) / 3, (2 * b003.z + b030.z) / 3);
564

    
565
            b201 = new Coordinate((2 * b300.x + b003.x - countScalarProduct(countDifferenceProduct(b003, b300), normalN1)
566
                                                         * normalN1.x) / 3, (2 * b300.y + b003.y - countScalarProduct(
567
                     countDifferenceProduct(b003, b300), normalN1)
568
                                                                                                   * normalN1.y) / 3,
569
                     (2 * b300.z + b003.z - countScalarProduct(countDifferenceProduct(b003, b300), normalN1) * normalN1.z) / 3);
570

    
571
            b102 = new Coordinate((2 * b003.x + b300.x - countScalarProduct(countDifferenceProduct(b300, b003), normalN3)
572
                                                         * normalN3.x) / 3, (2 * b003.y + b300.y - countScalarProduct(
573
                     countDifferenceProduct(b300, b003), normalN3)
574
                                                                                                   * normalN3.y) / 3,
575
                     (2 * b003.z + b300.z - countScalarProduct(countDifferenceProduct(b300, b003), normalN3) * normalN3.z) / 3);
576
            break;
577
         }
578
         case (6): {
579
            b210 = new Coordinate((2 * b300.x + b030.x) / 3, (2 * b300.y + b030.y) / 3, (2 * b300.z + b030.z) / 3);
580

    
581
            b120 = new Coordinate((2 * b030.x + b300.x) / 3, (2 * b030.y + b300.y) / 3, (2 * b030.z + b300.z) / 3);
582

    
583
            b021 = new Coordinate((2 * b030.x + b003.x) / 3, (2 * b030.y + b003.y) / 3, (2 * b030.z + b003.z) / 3);
584

    
585
            b012 = new Coordinate((2 * b003.x + b030.x) / 3, (2 * b003.y + b030.y) / 3, (2 * b003.z + b030.z) / 3);
586

    
587
            b201 = new Coordinate((2 * b300.x + b003.x) / 3, (2 * b300.y + b003.y) / 3, (2 * b300.z + b003.z) / 3);
588

    
589
            b102 = new Coordinate((2 * b003.x + b300.x) / 3, (2 * b003.y + b300.y) / 3, (2 * b003.z + b300.z) / 3);
590
         }
591
      }
592

    
593

    
594
      G = new Coordinate((b300.x + b030.x + b003.x) / 3, (b300.y + b030.y + b003.y) / 3, (b300.z + b030.z + b003.z) / 3);
595
      final Coordinate normalOfT = new Coordinate(0, 0, 1);
596

    
597
      A111 = new Coordinate((b300.x + b030.x + G.x) / 3, (b300.y + b030.y + G.y) / 3, (b300.z + b030.z + G.z) / 3);
598

    
599
      Coordinate e1 = normalizeVect(getNormal(1D / 2D, 0D));
600
      Coordinate e2 = normalizeVect(countDifferenceProduct(b120, b210));
601
      Coordinate normalOfPlane = setNormalVector(countCrossProduct(e1, e2), e2);
602

    
603
      A111 = new Coordinate(countProjectOnToPlane(b120, normalOfPlane, A111, normalOfT));//getNormal(4D/9D, 1D/9D)));
604

    
605
      B111 = new Coordinate((b003.x + b030.x + G.x) / 3, (b003.y + b030.y + G.y) / 3, (b003.z + b030.z + G.z) / 3);
606

    
607
      e1 = normalizeVect(getNormal(1D / 2D, 1D / 2D));
608
      e2 = normalizeVect(countDifferenceProduct(b012, b021));
609
      normalOfPlane = setNormalVector(countCrossProduct(e1, e2), e2);
610
      B111 = new Coordinate(countProjectOnToPlane(b012, normalOfPlane, B111, normalOfT));//getNormal(4D/9D, 4D/9D)));
611

    
612

    
613
      C111 = new Coordinate((b003.x + b300.x + G.x) / 3, (b003.y + b300.y + G.y) / 3, (b003.z + b300.z + G.z) / 3);
614

    
615
      e1 = normalizeVect(getNormal(0, 1 / 2D));
616
      e2 = normalizeVect(countDifferenceProduct(b201, b102));
617
      normalOfPlane = setNormalVector(countCrossProduct(e1, e2), e2);
618
      C111 = new Coordinate(countProjectOnToPlane(b201, normalOfPlane, C111, normalOfT));//getNormal(1D/9D, 4D/9D)));
619

    
620

    
621
      A201 = new Coordinate((b300.x + b210.x + b201.x) / 3, (b300.y + b210.y + b201.y) / 3, (b300.z + b210.z + b201.z) / 3);
622
      A021 = new Coordinate((b030.x + b120.x + b021.x) / 3, (b030.y + b120.y + b021.y) / 3, (b030.z + b120.z + b021.z) / 3);
623
      B102 = new Coordinate((b012.x + b102.x + b003.x) / 3, (b012.y + b102.y + b003.y) / 3, (b012.z + b102.z + b003.z) / 3);
624

    
625
      A102 = new Coordinate((A111.x + A201.x + C111.x) / 3, (A111.y + A201.y + C111.y) / 3, (A111.z + A201.z + C111.z) / 3);
626
      A012 = new Coordinate((B111.x + A021.x + A111.x) / 3, (B111.y + A021.y + A111.y) / 3, (B111.z + A021.z + A111.z) / 3);
627
      B201 = new Coordinate((C111.x + B111.x + B102.x) / 3, (C111.y + B111.y + B102.y) / 3, (C111.z + B111.z + B102.z) / 3);
628

    
629
      G = new Coordinate((A102.x + A012.x + B201.x) / 3, (A102.y + A012.y + B201.y) / 3, (A102.z + A012.z + B201.z) / 3);
630
   }
631

    
632

    
633
   /******************************************************************************************************************************
634
    * The method to print Bezier triangle to console
635
    */
636
   protected void printToConsole() {
637
      System.out.println("=========MACRO TRIANGLE=============================");
638
      System.out.println(b300.toString());
639
      System.out.println(b030.toString());
640
      System.out.println(b003.toString());
641
      System.out.println("Normals:");
642
      if (normalN1 != null) {
643
         System.out.println(normalN1.toString());
644
         System.out.println(normalN2.toString());
645
         System.out.println(normalN3.toString());
646
         System.out.println("Koeficients");
647

    
648
         System.out.println(b012.toString());
649
         System.out.println(b021.toString());
650
         System.out.println(b102.toString());
651
         System.out.println(b120.toString());
652
         System.out.println(b210.toString());
653
         System.out.println(b201.toString());
654
         //System.out.println(b111.toString());
655

    
656
      }
657
      System.out.println("======================================");
658
   }
659

    
660

    
661
   /******************************************************************************************************************************
662
    * The protected method for getting envelope of triangle
663
    * 
664
    * @return envelope of triangle
665
    */
666
   protected Envelope getEnvelope() {
667
      final Coordinate[] newPoint = new Coordinate[4];
668
      newPoint[0] = b003;
669
      newPoint[1] = b030;
670
      newPoint[2] = b300;
671
      newPoint[3] = b003;
672
      final CoordinateArraySequence newPointsTriangle = new CoordinateArraySequence(newPoint);
673

    
674
      final LinearRing trianglesPoints = new LinearRing(newPointsTriangle, new GeometryFactory());
675

    
676
      return trianglesPoints.getEnvelopeInternal();
677
   }
678

    
679

    
680
   /******************************************************************************************************************************
681
    * The protected method gets normal into Bezier triangle with barycentric coordinate
682
    * 
683
    * @param u -
684
    *                barycentric coordinate u
685
    * @param v -
686
    *                barycentric coordinate u
687
    * @return normal of surface
688
    */
689
   protected Bezier getBezierPatch(final int index) {
690
      switch (index) {
691
         case 0: {
692
            return new Bezier(b300, b030, G, b210, b120, A021, A012, A102, A201, A111);
693
         }
694
         case 1: {
695
            return new Bezier(b030, b003, G, b021, b012, B102, B201, A012, A021, B111);
696
         }
697
         case 2: {
698
            return new Bezier(b003, b300, G, b102, b201, A201, A102, B201, B102, C111);
699
         }
700
      }
701
      return null;
702
   }
703
}