low quality Bspline Curve get Using "points_to_bspline" function

[TsingQAQ]

Hi Paul,

When I'm using _points_to_bspline_ function to generate an airfoil curve, sometimes I get bad splines, especially at the trailing edge:

The problem only occurs when the following circumstances are both met:

1. when the function trys to generate a Bspline according to the data set of airfoil NACA 23012 ps: [naca 23012, ag 10, naca 4412] have been tested but only 23012 has the problem.
2. The y coordinates of the data set are not zero(explained bellow)

The data set are like follows: [NACA 23012 data sets] type: numpy array formate: [x_coordinates, y_coordinates, z_coordinates]

[[ 0.49889135 2. 0.08002602] [ 0.498087 2. 0.08022889] [ 0.4956829 2. 0.08083672] [ 0.49170192 2. 0.08183547] [ 0.4861903 2. 0.08320301] [ 0.47920573 2. 0.08491022] [ 0.47082381 2. 0.08692639] [ 0.46113132 2. 0.08921079] [ 0.45023623 2. 0.09172388] [ 0.43825518 2. 0.09442232] [ 0.42531681 2. 0.09726172] [ 0.41156117 2. 0.1001907 ] [ 0.39713804 2. 0.10316603] [ 0.38220515 2. 0.10613175] [ 0.36691723 2. 0.10903217] [ 0.35144693 2. 0.11181003] [ 0.33595523 2. 0.11441208] [ 0.32061128 2. 0.11677536] [ 0.3055793 2. 0.11884935] [ 0.29102024 2. 0.12058083] [ 0.27709016 2. 0.12192904] [ 0.26394262 2. 0.12285948] [ 0.25153484 2. 0.1233307 ] [ 0.2398575 2. 0.12313273] [ 0.22922185 2. 0.12210584] [ 0.21994481 2. 0.12025843] [ 0.21227763 2. 0.11774319] [ 0.20640394 2. 0.11480003] [ 0.20240754 2. 0.11169864] [ 0.20029478 2. 0.10864927] [ 0.2 2. 0.10579619] [ 0.20138579 2. 0.10328994] [ 0.20428678 2. 0.10120222] [ 0.20856322 2. 0.09941574] [ 0.21410069 2. 0.09779298] [ 0.22082608 2. 0.09618674] [ 0.22872648 2. 0.09448373] [ 0.23785485 2. 0.092601 ] [ 0.24830241 2. 0.09054866] [ 0.26016331 2. 0.08843892] [ 0.27320104 2. 0.08649423] [ 0.28709078 2. 0.08478819] [ 0.30167882 2. 0.08331917] [ 0.31680195 2. 0.08208039] [ 0.33229344 2. 0.08105932] [ 0.34798008 2. 0.080238 ] [ 0.36368539 2. 0.07959572] [ 0.37923863 2. 0.07910829] [ 0.39445812 2. 0.07876448] [ 0.40918481 2. 0.07853607] [ 0.42324626 2. 0.07841403] [ 0.4364929 2. 0.07837533] [ 0.44877347 2. 0.07841211] [ 0.45995633 2. 0.07849977] [ 0.46991115 2. 0.07862861] [ 0.47853369 2. 0.07877764] [ 0.48572 2. 0.07892887] [ 0.49139951 2. 0.07906737] [ 0.49549868 2. 0.07917852] [ 0.49797733 2. 0.07925075] [ 0.49880753 2. 0.07927448]]

I've done some tests to find some information about this problem:

• This fluctuation occurs at the direction perpendicular to the airfoil section plane.

• Not a display accuracy problem I'm using the follow command to get more accurat display, in adition I've exported to a step file to other commercial CADs（e.g CATIA）, still fine the problem. # increase the quality, as this shape has high curvature drawer = display.Context.DefaultDrawer().GetObject() drawer.SetDeviationAngle(drawer.DeviationAngle() / 100.) drawer.SetDeviationCoefficient(drawer.DeviationCoefficient() / 100.) thus I thought it was not a problem of display accuracy.

• May not be the NACA5 generating function problem I've tested two different ways to generate data sets of NACA 23012, though the result are different but both are not well: Using the builtin NACA5 airfoil generation fuction

Using the selig file(numpy.loadtxt) from profili(an Airfoil generation and analysis software)

• May be a problem of Numerical fluctuation Since the fluctuation direction is perpendicular to the airfoil plane, at first I thought it was a numerical accuracy problem(means that the numpy array at span location(which mean [:, 1] )like above may not be all equal), after I print those arrays it seems that they are all equal

Here's my way to generate an airfoil curve:

1. load data sets from Selig airfoil or from NACA generator function
2. reformat them to make a 3 dimension array: [[x_coordinate, y_coordinate:currently all 0, z_coordinate]...]
3. scale them due to the scale factor(only scale x and z coordinate)
4. offset them according to the leading edge point: (array in step3 + offset vector(a 3* 1 array)) [This step may be the main factor]
5. Using points_to_Bspline function to generate the curve

At step 4, if the y_coordinate of the offset vector is not zero, either int or float will cause such a problem(only NACA 23012 in my case), however, either x or z coordinate of the offset vector will not cause a similar problem.

I'm not sure whether this problem is caused by python or pythonocc.

[p-chambers]

Hi Will,

I am aware of this issue, however, I have only seen it if the points are scaled before doing the bspline interpolation. I am not sure that it issue is easily testable, unfortunately... even if I could write a test for it, there are too many airfoils in the Selig database for it to be appropriate to test all of them.

Have you tried using the standard primitives.Airfoil class constructor using the naca23012 from the Selig database (the bspline is transformed after the interpolation in this case)? Using the following code, I do not see the same issue with the curve:

if __name__ == '__main__':
    from airconics import primitives
    # Visualisation with Python-OCC (ensure plot windows are set to qt)
    from OCC.Display.SimpleGui import init_display
    display, start_display, add_menu, add_function_to_menu = init_display()

    LEPoint = [0., 2., 0.]
    ChordLength = 1
    Rotation = 0
    Twist = 0

    Af = primitives.Airfoil(LEPoint, ChordLength, Rotation, Twist,
                            SeligProfile='naca23012')

    display.DisplayShape(Af.Curve, update=True)
    display.DisplayShape(Af.ChordLine, update=True, color='black')
    start_display()

If not:

1. I've not implemented the bspline smoothing that was done in the Rhino version of airconics yet, so perhaps this could be a solution.

2. We might consider enforcing the curve to be 2D during interpolation, then convert to 3D before the transformation, as is the case in the pythonocc airfoil example.

Let me know if any of this works for you?

[TsingQAQ]

Hi Paul,

Some new test results:

• all the tested NACA 5 series have a similar problem [tested: 23012, 23009, 23015, 24012]

• I've also test the common research model CRM 65 airfoil, this time the problem is much the same like you said, if no scale the airfoil seems all right, but after scaling the curve is getting worse, I think smoothing may be a solution. (there are some unkown problem I cannot upload the figure CRM 65 curve)

About your suggestion:

1. Standard class: I've tried to use the standard builtin Airfoil class to get the naca 23012 Curve, no problem, though still, I want to find the reason of this problem,and since the Selig file is included in an egg file, I'm not able to get these files to see wether the origin data sets have any influnce on the final curve. Also, the NACA5 keyword of this function has not been implemented yet so I cannot test whether the airfoil function works all right.

2. Smoothing: I've awared that there's TODO mack about airfoil smoothing, I've no experience about the Rhino version of the code so I haven't try it.

3. 2D-3D: I've not tried this method yet, I think this may be a solution.

Actuatlly what really puzzled me is that even I'm using a 3D airfoil(numpy n* 3 array), after scaling, I can get a right naca23012 airfoil curve by doing an interpolating at this time if the y coordinates are all zero(means no span location), but onece I move them according to the span and then perform interpolating, the curve will not be good anymore, actually this problem could be avoided in mutiple ways like you suggested or interpolating first, but it is really strange...

[p-chambers]

Perhaps the backend bspline computation checks for if there is a y component in the points, and if not, it converts it to a 2d interpolation by default - I can't tell from the source code of GeomAPI_PointsToBspline or AppDef_BSplineCompute whether this is the case (I can't even find the source code for the latter, which I think does the computation). Or, it may be that the machine precision on the interpolation introduces some tiny differences in the y coordinate, and therefore it thinks the curve should extend in this plane to pass through those points.

In either case, the solution is to add a flag (or a separate function) for enforcing 2D interpolation for this case, and it would be more explicit to use that for generating the curve in the Airfoil class.

I'm a bit preoccupied preparing for a conference to do this any time soon. In the meanwhile, can I suggest that you create the BSpline using the Airfoil class, and allow the translation/scaling etc to be handled after the interpolation? Or you could have a go at writing the 2D function yourself (in which case, I would be happy to look at a pull request for this addition).

Sorry I can't be more helpful at the moment!

[TsingQAQ]

lolol Paul I think I may know the reason of this problem and need some time to proved that it is, just make sure points to bspline is interpolating, not fitting right?

[p-chambers]

I'm not sure what exact differences you are referring to, since I would say that interpolation and fitting roughly equate to the same thing in this case. It interpolates a cubic BSpline curve that approximately fits a set of points.

I provide a few flags to enable reuse for periodic curves, etc, but the default behaviour is to create a C^2 continuous cubic BSpline passing through the input points (within a tolerance), using GeomAPI_PointsToBspline.

[TsingQAQ]

Regret for my poor English, I mean 'smoothing' instead of fitting , interpolating means that the curve must go through all the given points, but smoothing means that the curve will not be forced to pass all the points if some points are very strange(e.g wrong points), it will use some algorithem(e.g least squares) to make the error small.

Some times smoothing curve will not tend to fluctuate like interpolating curve.(this may be a solution)

In this case,I find that if you set the degree of points to spline function to be 2, you can get a very good looking curve leaving anything unchanged. But if the degrees goes to 3, 4 everything is a mess, this may be a charactor of spline interpolating, I'm trying to figure out that.

[trelau]

Do the airfoils have blunt trailing edges or do they come to a sharp point? Are there sharp corners for the curve fitting routine to try and deal with? If so, it is common for curve interpolating and/or approximation to produce squiggly lines between points. At first glance this is what it looks like to me. A better way to fit an airfoil with a blunt TE is to first fit only the smooth parts of the airfoil and then close it out by adding an extra edge at the TE.

[TsingQAQ]

Yeah actually it has a blunt trailing edge and there's an extra line to close the trailing edge of the foil. Now my solution is just set the degree of bspline to be 2 to make a 'smooth looking' airfoil, but this shouldn't be a good solution since a degree 2 curve is not C2 continues.

[p-chambers]

Hi @TsingQAQ , Did you get to the bottom of this issue? I'm now back on this project, but can't seem to reproduce the problem: I tried to test with the NACA23012 points from airfoiltools, and then scaled the x and z and translated in the y, and the curve looks OK. I also tried with the xyz array you pasted above, and it still looks ok.

Could you provide a minimal working example?

[TsingQAQ]

Hi @p-chambers ,

Really glad to hear that you're back lolol. Still not, though I was not concentrating on this problem either. I've just changed the 'points_to_bspline' degree to be 2 to let it 'looks good' like said before and went over this.

Here I've test a DLR-F4 wing today, which looks like follows:

Here're some results: Set _points_tobspline max, and min degree = 2: Set _points_tobspline max, and min degree = 3: Set _points_tobspline max, and min degree = 4:

It clearly shows that the problem, though I've not went to the deep of the code thereafter and currently a bit hurry on a competition. The airfoil data process procedure(scale, translate... etc) has not been changed. Anything else I could provided for you?

This is the DLRF4 Section1 origin data: 202.13777 -17.56799 200.11639 -17.34976 198.09502 -17.17422 196.07364 -17.03 194.05224 -16.90811 192.03087 -16.81399 190.00949 -16.73368 187.98813 -16.667 185.96675 -16.61592 183.94537 -16.5844 181.92398 -16.56886 179.9026 -16.57147 177.88122 -16.5897 175.85985 -16.61834 173.83848 -16.66118 171.81711 -16.71618 169.79571 -16.78745 167.77434 -16.87738 165.75296 -16.97934 163.73158 -17.0944 161.7102 -17.22799 159.68884 -17.37848 157.66745 -17.549 155.64607 -17.74002 153.62469 -17.9511 151.60332 -18.17994 149.58194 -18.42497 147.56056 -18.6816 145.5392 -18.94064 143.51781 -19.20148 141.49643 -19.46119 139.47505 -19.7178 137.45367 -19.96813 135.4323 -20.2113 133.41093 -20.44651 131.38954 -20.67304 129.36816 -20.89351 127.34679 -21.10675 125.32541 -21.30808 123.30403 -21.50045 121.28265 -21.68313 119.26127 -21.8533 117.2399 -22.00953 115.21852 -22.15103 113.19714 -22.28131 111.17576 -22.40154 109.15439 -22.50523 107.13301 -22.60683 105.11163 -22.71053 103.09025 -22.76729 101.06888 -22.79315 99.0475 -22.79398 97.02612 -22.77386 95.00474 -22.75064 92.98337 -22.72141 90.96199 -22.68682 88.94061 -22.64643 86.91923 -22.59472 84.89786 -22.5347 82.87648 -22.46941 80.8551 -22.39609 78.83372 -22.31091 76.81235 -22.22135 74.79097 -22.12438 72.76959 -21.97722 70.74821 -21.81569 68.72682 -21.65283 66.70546 -21.48655 64.68408 -21.31631 62.6627 -21.13798 60.64133 -20.95422 58.61995 -20.76595 56.59856 -20.56644 54.57719 -20.35486 52.55582 -20.12804 50.53444 -19.89014 48.51306 -19.63895 46.49168 -19.36125 44.47031 -19.07124 42.44891 -18.76763 40.42755 -18.46143 38.40617 -18.15089 36.38479 -17.79624 34.36342 -17.42293 32.34204 -17.02933 30.32065 -16.61588 28.29927 -16.17728 26.27791 -15.70963 24.25653 -15.21803 22.23515 -14.69819 20.21377 -14.13207 18.19238 -13.5263 16.171 -12.85407 14.14963 -12.12773 12.12826 -11.34835 11.11757 -10.9203 10.10689 -10.46523 9.09619 -9.98374 8.59085 -9.73129 8.28764 -9.57507 8.08551 -9.46901 7.88336 -9.36123 7.68123 -9.25214 7.47909 -9.14161 7.27695 -9.0292 7.07481 -8.91469 6.87268 -8.79816 6.67053 -8.67954 6.4684 -8.55895 6.26626 -8.43603 6.06413 -8.31028 5.86198 -8.1821 5.65985 -8.05215 5.45772 -7.92106 5.25557 -7.78838 5.05343 -7.65157 4.8513 -7.51161 4.64915 -7.36953 4.44702 -7.22589 4.24489 -7.07865 4.04275 -6.92586 3.84061 -6.76732 3.63847 -6.6025 3.43634 -6.43088 3.23419 -6.25189 3.03206 -6.0653 2.82992 -5.86726 2.62778 -5.65546 2.42564 -5.43273 2.22351 -5.20203 2.02138 -4.96658 1.81923 -4.7289 1.61709 -4.48296 1.41496 -4.22796 1.21281 -3.95506 1.01068 -3.65506 0.80855 -3.3175 0.6064 -2.9275 0.40426 -2.46056 0.20213 -1.84056 0.10106 -1.35056 0.05052 -1.01056 0.0202 -0.74056 0.0101 -0.47056 0 0 0.0101 0.68039 0.0202 0.98945 0.05052 1.35329 0.10106 1.78288 0.20213 2.28216 0.40426 3.02196 0.6064 3.53135 0.80855 3.95182 1.01068 4.31566 1.21281 4.63502 1.41496 4.91803 1.61709 5.16665 1.81923 5.38294 2.02138 5.57498 2.22351 5.74275 2.42564 5.90532 2.62778 6.05895 2.82992 6.20145 3.03206 6.33487 3.23419 6.4612 3.43634 6.58047 3.63847 6.69366 3.84061 6.80079 4.04275 6.90388 4.24489 6.9999 4.44702 7.09086 4.64915 7.17931 4.8513 7.26299 5.05343 7.34215 5.25557 7.41718 5.45772 7.48832 5.65985 7.5559 5.86198 7.62147 6.06413 7.68627 6.26626 7.74988 6.4684 7.81177 6.67053 7.87121 6.87268 7.92985 7.07481 7.98795 7.27695 8.04522 7.47909 8.10117 7.68123 8.15519 7.88336 8.20739 8.08551 8.2582 8.28764 8.30772 8.59085 8.37957 9.09619 8.49291 10.10689 8.70557 11.11757 8.90179 12.12826 9.08082 14.14963 9.38099 16.171 9.63244 18.19238 9.83091 20.21377 10.00045 22.23515 10.14928 24.25653 10.26461 26.27791 10.35216 28.29927 10.41979 30.32065 10.46984 32.34204 10.49737 34.36342 10.50461 36.38479 10.49836 38.40617 10.47362 40.42755 10.43543 42.44891 10.38389 44.47031 10.31928 46.49168 10.2394 48.51306 10.15011 50.53444 10.04912 52.55582 9.93592 54.57719 9.81458 56.59856 9.68435 58.61995 9.5458 60.64133 9.39834 62.6627 9.2405 64.68408 9.07886 66.70546 8.91112 68.72682 8.73531 70.74821 8.55237 72.76959 8.36234 74.79097 8.16518 76.81235 7.96316 78.83372 7.75281 80.8551 7.5366 82.87648 7.31513 84.89786 7.0868 86.91923 6.85338 88.94061 6.61372 90.96199 6.36915 92.98337 6.11862 95.00474 5.86387 97.02612 5.60313 99.0475 5.33825 101.06888 5.06843 103.09025 4.79361 105.11163 4.51466 107.13301 4.23069 109.15439 3.94149 111.17576 3.64734 113.19714 3.34828 115.21852 3.04509 117.2399 2.7358 119.26127 2.42156 121.28265 2.10218 123.30403 1.77759 125.32541 1.44808 127.34679 1.1126 129.36816 0.76991 131.38954 0.42118 133.41093 0.06427 135.4323 -0.2996 137.45367 -0.67123 139.47505 -1.04995 141.49643 -1.43582 143.51781 -1.82962 145.5392 -2.23236 147.56056 -2.64282 149.58194 -3.06216 151.60332 -3.48529 153.62469 -3.91287 155.64607 -4.3501 157.66745 -4.79476 159.68884 -5.24579 161.7102 -5.70622 163.73158 -6.17156 165.75296 -6.64565 167.77434 -7.14817 169.79571 -7.62269 171.81711 -8.11445 173.83848 -8.62238 175.85985 -9.13387 177.88122 -9.66474 179.9026 -10.20621 181.92398 -10.75218 183.94537 -11.32472 185.96675 -11.8808 187.98813 -12.46659 190.00949 -13.05736 192.03087 -13.66121 194.05224 -14.25943 196.07364 -14.89217 198.09502 -15.52844 200.11639 -16.17024 202.13777 -16.82185

[p-chambers]

@TsingQAQ , I recreated the issue you mention. The curvature does not extend into the y dimension, so contrary to what I said previously, I don't think enforcing a 2D interpolation will fix the issue.

However, we could use some of the internal workings of GeomAPI_PointsToBSpline to enforce a smooth airfoil curve. I messed around with the GeomAPI_PointsToBSpline method that uses a weighted variational smoothing algorithm, but this required an unnecessarily high curve degree to successfully produce the shape. I had more luck using the optional tangents argument of points_to_bspline, enforcing a tangent vector on the start and end of the curve, using the first and last two points from the airfoil data. Here's the code:

if __name__ == '__main__':
    from airconics import primitives
    import numpy as np
    import airconics.AirCONICStools as act
    # Visualisation with Python-OCC (ensure plot windows are set to qt)
    from OCC.Display.SimpleGui import init_display
    display, start_display, add_menu, add_function_to_menu = init_display()
    drawer = display.Context.DefaultDrawer().GetObject()
    drawer.SetDeviationAngle(drawer.DeviationAngle() / 100.)
    drawer.SetDeviationCoefficient(drawer.DeviationCoefficient() / 100.)

    LEPoint = [0., 1., 0.]
    Af = primitives.Airfoil(LEPoint)

    # I copy and pasted your points above into a file called DLR_F4.txt...
    points = np.loadtxt('DLR_F4.txt')
    # Give the points a non zero y, in case this is causing the problem
    points = np.column_stack([points[:, 0], np.ones_like(points[:, 0]),
                              points[:, 1]])

    # Create the curve without tangents, showing poor Trailing Edge:
    curve1 = act.points_to_bspline(points, deg=3).GetObject()
    display.DisplayShape(curve1, color='blue')
    # plot the control points, showing why the curve is poor at the TE
    for i in range(curve1.NbPoles()):
        display.DisplayShape(curve1.Pole(i), color='black', update=True)

    # Create the red curve enforcing the upper and lower trailing edge curves
    # to have tangent vectors in the direction of the adjacent points
    tangents = np.array([points[1, :] - points[0, :],
                        points[-1, :] - points[-2, :]])

    curve2 = act.points_to_bspline(points, deg=3, tangents=tangents)

    display.DisplayShape(curve2, update=True, color='red')

    start_display()

Here's a screen grab of the trailing edge:

Does this solve your problem? I don't think that this is an error with points_to_bspline since it is not specifically created for producing airfoils, but this may a useful addition to the Airfoil class.

[TsingQAQ]

Hi @p-chambers,

many thanks for your detailed advice sorry for the late reply I was just had enough time to test this. Giving a tangent constraint should work well. Though, there're some changes on the return of the _points_tobspline and may cause new problem, I've tried to fix it but could not work it out.

just at the end of your code. I'm using a make edge function to close the trailing edge like bellow: which should looks like this if we do not use tangent constraint kwargs: The blue one is the trailing edge, it will close the section.

If we give a tangent constarint at the _points_tobspline, function the trailing edge will be the following:

They all use the same codes bellow except for the tangent kwargs:

    # your code above....
    curve2 = act.points_to_bspline(points, deg=3, tangents=tangents)

    # close the trailing edge:
    crv = curve2.GetObject()
    te_edge = act.make_edge(crv.Value(1), crv.Value(0))

    display.DisplayShape(curve2, update=True, color='red')
    display.DisplayShape(te_edge, update=True, color='blue')

    start_display()

Seems that the return curve have changed, the start and end points(use .Value(0) and .Value(1) to specify this) of the return curve no longer locates at the trailing edges, so I wonder how could I close the trailing edge if using a tangent constraint.

[p-chambers]

If you replace your te_edge line to use the StartPoint and EndPoint methods, you should get the desired effect:

    te_edge = act.make_edge(crv.StartPoint(), crv.EndPoint())

The coordinate system along the curve seems to behave differently when specifying the tangent vector, in that the maximum Value of the curve is no longer 1 (you can see this with crv.LastParameter()). This may have something to do with subtle differences between the underlying OCC interpolation methods used when tangents are given (GeomAPI_PointsToBSpline vs GeomAPI_Interpolate).

I don't understand why this is the case with OCC, but at least we know what the issue is.

[TsingQAQ]

Wow It works! Really thanks Paul finally looks like a plan.

Shall we close this issue then?

[p-chambers]

Excellent. I'd like this to be an optional method in the Airfoil class, so I'll raise a new issue for that and close this one :).