[SOLVED] Interpolation

[mari_ia]

Hi there!

Could you please help me to resolve the error appeared as "External: FLT INVALID OPERATION" (on the line: Sarray[n,1]:= Iarray[m,1]+((Iarray[m+1,1]-Iarray[m,1])/(Iarray[m+1,0]-Iarray[m,0]))*(Sarray[n,0]-Sarray[n-1,0])?

My task is to do the interpolation of experimental data. There are X and Y points that I put in two columns of StringGrid1 table. Then I transfer these data to the array inter1. After that, I divide the segments between two neighboring points in half and calculate the value of Y corresponding to this point. In such a way, I fill in the array inter2. There may be several cycles of dividing the segments in half. At the final stage, I transfer the elements of inter2 to StringGrid2.

Seems I do not take into account the change in the rows count of the array inter2. Maybe this is the reason.

Please find the fragment of my code below. I will be really glad for your comments!

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

1. type
2.   InputArr = array of array of real;
3.   SplitArr = array of array of real;
4.   OutputArr = array of array of real;
5.  
6.   { TForm1 }
7.  
8.   TForm1 = class(TForm)
9.     // ...
10.  
11.   private
12.  
13.   public
14.   Ni, rsplit : integer;
15.   inter1 : array of array of real;
16.   inttemp : array of array of real;
17.   inter2 : array of array of real;
18.   procedure Splitting(var Iarray : InputArr; Sarray : SplitArr; rinput : integer);
19.   procedure Outputs(Arr : OutputArr; ncol : integer);
20.   end;
21.  
22. var
23.   Form1: TForm1;
24.  
25. implementation
26.  
27. {$R *.lfm}
28.  
29. { TForm1 }
30.  
31. {0a - number of input rows}
32. procedure TForm1.Edit1Exit(Sender: TObject);
33. begin
34.   Ni:=StrToInt(Edit1.Text);
35.   StringGrid1.RowCount:=Ni+1;
36. end;
37.  
38. procedure TForm1.Splitting(var Iarray : InputArr; Sarray : SplitArr; rinput : integer);
39. var m, n : integer;
40. begin
41.   rsplit:=2*rinput-1;
42.   Edit2.Text:=IntToStr(rsplit);
43.   //SetLength(Sarray, rsplit, 2);
44.   m:=0;
45.   while m <= rinput-1 do
46.   begin
47.    n:=2*m;
48.    Sarray[n,0]:=Iarray[m,0];
49.    Sarray[n,1]:=Iarray[m,1];
50.    if m = rinput-1 then
51.    m:=m+1
52.    else
53.    begin
54.    n:=n+1;
55.    Sarray[n,0]:=(Iarray[m,0]+Iarray[m+1,0])/2;
56.    Sarray[n,1]:= Iarray[m,1]+((Iarray[m+1,1]-Iarray[m,1])/(Iarray[m+1,0]-Iarray[m,0]))*(Sarray[n,0]-Sarray[n-1,0]);
57.    m:=m+1;
58.    end;
59.    end;
60. end;
61.  
62. procedure TForm1.Outputs(Arr : OutputArr; ncol : integer);
63. var rout, cout : integer;
64. begin
65.    //Processed Array -> StringGrid2
66.   for cout:=0 to 1 do
67.    begin
68.     for rout:=0 to StringGrid2.RowCount-2 do
69.     begin
70.     if cout = 0 then StringGrid2.Cells[cout,rout+1]:=FloatToStr(Arr[rout,cout])
71.     else StringGrid2.Cells[ncol,rout+1]:=FloatToStr(Arr[rout,cout]);
72.     end;
73.    end;
74. end;
75.  
76. {1 - interpolation}
77. procedure TForm1.Button1Click(Sender: TObject);
78. var
79.   i, j, k, ncycle, acol, qnt, No : integer;
80.  
81. begin
82.   Chart1LineSeries1.Clear;
83.   Chart1LineSeries2.Clear;
84.  
85.   //No:=StrToInt(Edit2.Text);
86.   //StringGrid2.RowCount:=No+1;
87.  
88.   SetLength(inter1, StringGrid1.RowCount-1, 2);
89.   //SetLength(inter2, StringGrid2.RowCount-1, 2);
90.   acol:=1;
91.  
92.   //checking the elements
93.   for k:=1 to StringGrid1.RowCount-1 do
94.     begin
95.      if (StringGrid1.Cells[0, k]='') or (StringGrid1.Cells[1, k]='')
96.      then
97.       begin
98.      ShowMessage('Check your input data');
99.      Exit;
100.       end;
101.     end;
102.  
103.   //StringGrid1 -> Input Array
104.   for i:=0 to StringGrid1.RowCount-2 do
105.    begin
106.     for j:=0 to 1 do
107.     begin
108.     inter1[i,j]:= StrToFloat(StringGrid1.Cells[j,i+1]);
109.     end;
110.    end;
111.  
112.   ncycle:=StrToInt(Edit2.Text) Div StrToInt(Edit1.Text);
113.   if ncycle > 0 then
114.      if ncycle*StrToInt(Edit1.Text)/StrToInt(Edit2.Text) <= 0.80
115.      then ncycle:=ncycle+1 else ncycle:=ncycle+0;
116.  
117.   if ncycle <= 0 then begin
118.   ShowMessage('Check the number of output data points');
119.   Exit;
120.   end;
121.  
122.   //Input Array -> Interpolated Array
123.  
124.   for qnt:= 0 to ncycle-1 do
125.   begin
126.    if qnt = 0 then
127.     begin
128.     SetLength(inter2, 2*Ni-1, 2);
129.     Splitting(inter1, inter2, Ni);
130.     end
131.    else
132.     begin
133.     SetLength(inter2, rsplit, 2);
134.     Splitting(inter2, inter2, rsplit);
135.     end;
136.   end;
137.   StringGrid2.RowCount:=rsplit+1;
138.  
139.   //Interpolated Array -> SringGrid2
140.   Outputs(inter2, acol);
141.  
142. end;

[wp]

At first you should help us and format your code for better readability. Being the author of your post you should see a "Modify" button above your post. Click on it to get your text back into editing mode. Then select, with the mouse, the code, beginning with "type" and ending with "end;". Above the editor you'll find several combo boxes: click on the "Code" combobox and select "pascal". This surrounds your code block with [ code ] tags. Check with the "Preview" button if everything is correct, and finally click "Post".

Since this is quite long code I'd also mention that usually it is better if you submit a full compilable project. Copy your project to a new folder, and in the copy remove everything not needed. Check to make sure that the issue still exists in the minimized project. Then pack the .pas, .lfm, .lpi and .lpr files (as well as any data files if needed) into a shared zip which you can upload under "Attachments and other options". Do not include the exe, ppu or other compiler-generated files because there is an upload limit of only 500 KB, IIRC.

[avra]

Interpolation is easy. You are probably over complicating things. You can calculate in runtime any f(x) point between neighbour points f(x1) and f(x2). Do you really need cutting in half each time and creating new array of points?

p.s. I did not look into your code, just your text


This might be interesting to you: https://bitbucket.org/avra/interpolate1234d/src/master/

[Thaddy]

Why not simply explain it like the sum of the values divided by the number of values?
This is not even college grade but basic from 5th grade in the Netherlands ( so from 9-10 years old)
So, given an array [0.2] and values being 1,2,3 you divide the sum  of 1+2+3 by  3.
 
Avra is right but he uses a notation that juniors are not familiar with and I suspect you are a junior.

More complex interpolation, for example for image processing, use a 3 dimension plain, not a 2 dimension plain as I described here

[Curt Carpenter]

Code is a little hard to follow.  But if you start with n points (and I did the math right), after m interations you should have (2^m)*(n-1)+1 points.  Maybe that can help you find out what's not working.  Arrays are growing fast at any rate.

[avra]

More complex interpolation, for example for image processing, use a 3 dimension plain, not a 2 dimension plain as I described here

In Interpolate1234D I cover 1D, 2D, 3D and 4D. I didn't have a need for more dimensions 

[mari_ia]

Hi:)

I did:

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

1. function TForm1.Splitting(Iarray : InputArr; rinput : integer) : SplitArr;
2. var m, n : integer;
3. begin
4.   rsplit:=2*rinput-1;
5.   Edit2.Text:=IntToStr(rsplit);
6.   SetLength(Result, rsplit, 2);
7.   m:=0;
8.   while m <= rinput-1 do
9.   begin
10.    n:=2*m;
11.    Result[n,0]:=Iarray[m,0];
12.    Result[n,1]:=Iarray[m,1];
13.    if m = rinput-1 then
14.    m:=m+1
15.    else
16.    begin
17.    n:=n+1;
18.    Result[n,0]:=(Iarray[m,0]+Iarray[m+1,0])/2;
19.    Result[n,1]:= Iarray[m,1]+((Iarray[m+1,1]-Iarray[m,1])/(Iarray[m+1,0]-Iarray[m,0]))*(Result[n,0]-Result[n-1,0]);
20.    m:=m+1;
21.    end;
22.    end;
23. end;
24.  
25. procedure TForm1.Button1Click(Sender: TObject);
26. //...
27. //Input Array -> Interpolated Array
28.   for qnt:= 0 to ncycle-1 do
29.   begin
30.    if qnt = 0 then
31.     begin
32.     SetLength(inter2, 2*Ni-1, 2);
33.     inter2:=Splitting(inter1, Ni);
34.     end
35.    else
36.     begin
37.     SetLength(inttemp, 2*rsplit-1, 2);
38.     inttemp:=Splitting(inter2, rsplit);
39.     SetLength(inter2, rsplit, 2);
40.     inter2:=inttemp;
41.     end;
42.   end;
43.   StringGrid2.RowCount:=rsplit+1;
44. //...
45.  

[Paolo]

Hello,

if you known a priopri the number of cycles you can dimension the oputput array size at the begin, avoding continuosly resizing the array length. (I assuming linear interpolation amgon cardinal points).

Here my proposal (not in deep tested) and maybe not optimal one...

It works even if the starting sampling is not with constant Dx. You can easily adpat the code at your case (input section data)

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

1. procedure TForm1.Button2Click(Sender: TObject);
2. var
3.   NPnt, NCycls : integer;
4.   Ofs, k, i, j : integer;
5.   DX, DY, M, DXk : double;
6. begin
7.   //--- INPUT dimension
8.   NPnt:=6;
9.   NCycls:=4;
10.   //--- DIMENSIONING OUTPUT ARRAY
11.   k:=trunc(Power(2, NCycls));
12.   SetLength(ArrayXY, (NPnt-1)*k+1);
13.   //--- set cardinal points (X,Y) from input values
14.   for i:=0 to NPnt-1 do begin
15.     ArrayXY[i*k].X:=i;
16.     ArrayXY[i*k].Y:=2*i;
17.   end;
18.   //--- fill sub-set inside any sub-interval with linear interpolation
19.   for i:=0 to NPnt-2 do begin
20.     Ofs:=i*k;                              //current sub-set start index
21.     DX:=ArrayXY[Ofs+k].X-ArrayXY[Ofs].X;   //current sub-set DX
22.     DY:=ArrayXY[Ofs+k].Y-ArrayXY[Ofs].Y;   //current sub-set DY
23.     M:=DY/DX;
24.     //--- loop on points inside the sub-interval (extreme points are the cardinal points, already available in the array)
25.     DXk:=DX/k;                             //inter-sub-interval Delta-X
26.     for j:=1 to k-1 do begin
27.       ArrayXY[Ofs+j].X:=ArrayXY[Ofs].X+j*DXk; //j-esimo element in sub-interval i-esimo
28.       ArrayXY[Ofs+j].Y:=(ArrayXY[Ofs+j].X-ArrayXY[Ofs].X)*M+ArrayXY[Ofs].Y; //linear interpolation
29.     end;
30.   end;
31.   Listbox1.clear;
32.   for i:=0 to Length(ArrayXY)-1 do
33.     ListBox1.Items.Add('X = '+ArrayXY[i].X.Tostring+', Y = '+ArrayXY[i].Y.Tostring);
34. end;
35.