Questions about TFuncSeries

[hedgehog]

Hi!

I want to improve the accuracy of the results of plotting a function containing sharp peaks.

For example:

Code: Pascal  [Select][+][-]

1. procedure TForm1.Chart1FuncSeries1Calculate(const AX: Double; out AY: Double);
2. begin
3.   AY:= 1/(abs(AX)+0.001);
4. end;

Obviously, the peak height at 0 should be 1000.
However, in the attached screenshot, you can see that the peak height is less than 500 (with step=1).

To do this, I need to change the code of the procedure TDrawFuncHelper.ForEachPoint

Question: What's the best way to do this without modifying the TAChart package?

Create your own series? (TMyFuncSeries)

[Thaddy]

It is not that the peak isn't ~1000, it is because the scaling is such that it won't draw it because not enough samples are in the display range on X axis. You have to adapt the drawing such, that given the scaling and the screen resolution, a peak is drawn at ~1000 regardless:
For the purpose of display, e.g anything above ~500 - ~800 (dpe3nds on steepness) should be scaled to a higher value or simply to the peak straight away.
That way a peak value is always displayed.
You can also use an envelope follower to smooth out the display range vs the real sample range. (Very easy to implement: a running queue (over X-axis) with a running average of samplerate/displayX elements + a bias > 1.

[wp]

In this special case, you could check whether AX is close enough to zero. "Close enough" means: The value 0 is in the interval which is spanned on the x axis by one TFuncSeries.Step. You must calculate the size of Step on the x axis and then check the agreement by means of the functions IsZero or SameValue from the Math unit where this step size can be used as tolerance variable.

Code: Pascal  [Select][+][-]

1. // Calculate the size of Step in Graph units
2. // Assuming here that graph units are the same as axis units, i.e. there is
3. // no transformation involved.
4. function TForm1.StepToGraph(X: Double; AStep: Integer): Double;
5. var
6.   imgX: Integer;
7.   nextX: Double;
8. begin
9.   imgX := Chart1.XGraphToImage(X);
10.   nextX := Chart1.XImageToGraph(imgX + AStep);    // the value after a "step"
11.   Result := nextX - X;
12. end;
13.  
14. // Save the result of StepToGraph in some variable, e. g., FGraphStep
15. procedure TForm1.Chart1FuncSeries1Calculate(const AX: Double; out AY: Double);
16. begin
17.   if IsNaN(FGraphStep) then
18.     FGraphStep := StepToGraph(0.0, Chart1FuncSeries1.Step);
19.   if IsZero(AX, FGraphStep) then
20.     AY := 1/0.001
21.   else
22.     AY:= 1/(abs(AX)+0.001);
23. end;
24.  
25. // Enforce recalculation of FGraphStep when the chart is resized.
26. procedure TForm1.Chart1Resize(Sender: TObject);
27. begin
28.   FGraphStep := Nan;
29. end;
30.  
31. // Enforce recalculation of FGraphStep when the chart is zoomed.
32. procedure TForm1.Chart1ExtentChanging(ASender: TChart);
33. begin
34.   FGraphStep := NaN;
35. end;
36.  
37. procedure TForm1.FormCreate(Sender: TObject);
38. begin
39.   FGraphStep := NaN;
40. end;        [/quote]
41. This code calculates the size of one "step" in the x axis "graph" units and stores the result in variable FGraphStep. It triggers recalculation of FGraphStep when the scaling of the plot changes, i.e. upon resizing and zooming.

[Thaddy]

Ignore my answer: as per wp's demo there is already code to do similar.

[hedgehog]

Hi.

WP: Your answer is good, but it only works in this specific case!
For any other function, we won't get the correct graph.

Actually, I was trying to implement a polyline simplification algorithm for TChart/TFuncSeries.
But after two hours of studying TChart's internal mechanisms, I realized it was quite complicated.

[wp]

I realized it was quite complicated.

Yes, it is...

Your answer is good, but it only works in this specific case!
For any other function, we won't get the correct graph.

I'd say there's a infinite number of math functions, it is impossible to take care of all difficulties in a general solution. Can you explain what is common to the functions that you are trying to plot?

a polyline simplification algorithm

You mean you want to create an array of points to be drawn by Canvas.PolyLine? This would have the advantage that, after creation of the array, you can insert the missing data points at the exact location of the "pseudo singularities". The disadvantage is that you must recalculate the array explicitely whenever the chart is zoomed or enlarged. You could even use a plain old TLineSeries for this...

[hedgehog]

I'd say there's a infinite number of math functions, it is impossible to take care of all difficulties in a general solution. Can you explain what is common to the functions that you are trying to plot?

What all these functions have in common is that they may have sharp peaks. Or they may not.

You mean you want to create an array of points to be drawn by Canvas.PolyLine?

Yes.
In the attached project, I attempted to create a universal procedure for plotting a function graph.

This should always work (at least, I think so).

[wp]

Interesting. I must confess that I do not understand this algorithm, without having looked into it too deeply... Can you give a verbal description?

Certainly, a patch for integration into TFuncSeries would be welcome. But maybe this method should be available as a separate option because after adding a counter into your calc function I see that there are almost 63,000 function calls to draw the curve in your demo - which may be critical when calculation of the curve to be plotted is complicated. Two mutually exclusive features: current TFuncSeries: many data points, a "few" function calls, vs. your method: very few data points, but very many function calls.

A requirement for adding to TFuncSeries would be also that it must be compatible with the DomainExclusion feature of TFuncSeries, i.e. it must respect the regions in which the function is not defined (e.g. sqrt(x) with domain exclusion for x < 0).

[hedgehog]

Oh, yes, my apologies. The "swing door trend" algorithm is used here.

It's a very simple and undemanding algorithm for simplifying polylines.
Surprisingly, despite its widespread use, it doesn't have a Wikipedia entry.

Here is a very simple description of the algorithm:

https://softwaredocs.weatherford.com/cygnet/94/Content/Topics/History/CygNet%20Swinging%20Door%20Compression.htm

But maybe this method should be available as a separate option because after adding a counter into your calc function I see that there are almost 63,000 function calls to draw the curve in your demo - which may be critical when calculation of the curve to be plotted is complicated.

The algorithm is very simple. Each step performs several multiplications of real numbers and several comparisons.
For any modern processor this is not a problem.

A requirement for adding to TFuncSeries would be also that it must be compatible with the DomainExclusion feature of TFuncSeries, i.e. it must respect the regions in which the function is not defined (e.g. sqrt(x) with domain exclusion for x < 0).

Oh, well, this one is simple. You need to plot the function graph in chunks, not the entire thing.

[wp]

The "swing door trend" algorithm is used here.

It's a very simple and undemanding algorithm for simplifying polylines.
Surprisingly, despite its widespread use, it doesn't have a Wikipedia entry.

Here is a very simple description of the algorithm:

https://softwaredocs.weatherford.com/cygnet/94/Content/Topics/History/CygNet%20Swinging%20Door%20Compression.htm

I'll take a closer look.

But maybe this method should be available as a separate option because after adding a counter into your calc function I see that there are almost 63,000 function calls to draw the curve in your demo - which may be critical when calculation of the curve to be plotted is complicated.

The algorithm is very simple. Each step performs several multiplications of real numbers and several comparisons.

No, I'm not talking about the algorithm itself, but about the calls to MyFunc. In your case, it is simple, but I could imagine that such a function can be quite complex. To count the function calls you can introduce some counting variable which you increment in MyFunc, reset at the top of Paintbox OnPaint handler and read at the end of it:

Code: Pascal  [Select][+][-]

1. var
2.   NumCalls: Integer;
3.  
4. procedure TForm1.PaintBox1Paint(Sender: TObject);
5. begin
6.   NumCalls := 0;
7.   // ... your code here ...
8.   Caption := IntToStr(NumCalls);
9. end;
10.  
11. function TForm1.MyFunc(x: double): double;
12. begin
13.   Result:= -1/(abs(x)+0.001);
14.   inc(NumCalls);
15. end;

                           

[hedgehog]

No, I'm not talking about the algorithm itself, but about the calls to MyFunc. In your case, it is simple, but I could imagine that such a function can be quite complex.

Just need to lower the sampling frequency.
Instead:

Code: Pascal  [Select][+][-]

1. dx:= kx/100; //100 samples/pixel

need:

Code: Pascal  [Select][+][-]

1. dx:= kx/15; //15 samples/pixel

[wp]

Are you sure that your left-right vector math is correct? Assume vector A pointing to the right, A = (1, 0), and B pointing upward B = (0, 1).

Code: [Select]

   ^ B
   |
   |
     - - -  > A
Looking along the direction of A the left vector is B, and conversely, the same result is obtained when looking along B where A is on the right side, i.d. the "left" vector must be B.

Doing your calculation in GetLeftVector, however, returns A as the left vector:

Code: Pascal  [Select][+][-]

1. program project1;
2.  
3. {$modeswitch advancedrecords}
4.  
5. type
6.   TVector = record
7.     x, y: double;
8.     function IsLeft(V: TVector): boolean;
9.   end;
10.  
11. function TVector.IsLeft(V: TVector): boolean;
12. begin
13.   Result:= (x*V.y) > (V.x*y);
14. end;
15.  
16. function Vector(x, y: double): TVector; inline;
17. begin
18.   Result.x:= x;
19.   Result.y:= y;
20. end;
21.  
22. function GetLeftVector(A, B: TVector): TVector; inline;
23. begin
24.   if A.IsLeft(B) then Result:= A else Result:= B;
25. end;
26.  
27. procedure WriteVector(AName: String; v: TVector);
28. begin
29.   WriteLn(AName, ' = (', v.x:0:3, ',', v.y:0:3, ')');
30. end;
31.  
32. var
33.   A, B: TVector;
34. begin
35.   A := Vector(1, 0);
36.   B := Vector (0, 1);
37.  
38.   WriteVector('A', A);
39.   WriteVector('B', B);
40.   WriteVector('GetLeftVector(A, B)', GetLeftVector(A, B));
41.  
42.   ReadLn;
43. end.

[wp]

No, I'm not talking about the algorithm itself, but about the calls to MyFunc. In your case, it is simple, but I could imagine that such a function can be quite complex.

Just need to lower the sampling frequency.
Instead:

Code: Pascal  [Select][+][-]

1. dx:= kx/100; //100 samples/pixel

need:

Code: Pascal  [Select][+][-]

1. dx:= kx/15; //15 samples/pixel

Ah, as I understand now you are doing some massive oversampling of each pixel which helps you to detect the sharp maximum, and use the "swing door trend" algorithm to reject all duplicate pixels and pixels which do not change in y significantly.

Although the peak detection is greatly improved compared with the current TAChart routine calculating only a single value per pixel it is not "fool-proof": If I am stupid enough to plot your function y = 1/(abs(x)+0.001) in the range -1E9 to +1E9 there are also cases, where the position of the maximum is missed visibly.

[hedgehog]

Ah, as I understand now you are doing some massive oversampling of each pixel which helps you to detect the sharp maximum, and use the "swing door trend" algorithm to reject all duplicate pixels and pixels which do not change in y significantly.

Not exactly. It's not about finding a maximum. The SDT algorithm replaces a group of points with a line, guaranteeing that all "discarded" points are no more than epsilon away from the line.

I tried plotting it over the range of -1e9 to 1e9.
The result, as you can see in the image, isn't too bad.
The peak height is certainly smaller, but it's still easily discernible visually.

[wp]

Ah, as I understand now you are doing some massive oversampling of each pixel which helps you to detect the sharp maximum, and use the "swing door trend" algorithm to reject all duplicate pixels and pixels which do not change in y significantly.

Not exactly. It's not about finding a maximum. The SDT algorithm replaces a group of points with a line, guaranteeing that all "discarded" points are no more than epsilon away from the line.

The algorithm, ok. But applied to the topic of this post, massive oversampling is the only way to detect such a sharp maximum (or minimum), isn't it?

I tried plotting it over the range of -1e9 to 1e9.
The result, as you can see in the image, isn't too bad.
The peak height is certainly smaller, but it's still easily discernible visually.

I also client-aligned the Paintbox and increased/decreased the width of the form. Then the peak height was jumping up and down.

One more question:
Wouldn't it be sufficient to compare the slopes of the test lines, rather than doing these multiple left-right comparisons? Of course, division is more expensive than multiplication. But how many multiplications have the same effect as one division?