Expiring Today—Celebrate National IT Professionals Day with 3 months of free Premium Membership. Use Code ITDAY17

x
Solved

# Interpolation, approximation?

Posted on 2002-04-04
Medium Priority
378 Views
Hi all.

Let say I have three base points (P1, P2, P3) on surface. Every base point defined by X1, X2, X3 and
Y1, Y2, Y3 coordinates. Every base point has one associated value Vn (V1, V2, V3). I have to calculate
approximated value V4 for point P4 defined by coordinates X4 and Y4. It's just a simplification of the

In deed, I need a fast algorithm that allow me to calculate value for any point on the surface where
number of base points may be different, let say N.

type
TSFPoint = record
X: Integer;
Y: Integer;
V: Integer;
end;

TSFPoints = array of TSFPoint;

implementation

function GetValueFor(X,Y: Integer; Points: TSFPoints): Integer;
begin
?????
end;
-----
Igor.

0
Question by:ITugay
[X]
###### Welcome to Experts Exchange

Add your voice to the tech community where 5M+ people just like you are talking about what matters.

• Help others & share knowledge
• Earn cash & points
• Learn & ask questions
• 9
• 6
• 3

LVL 1

Expert Comment

ID: 6920199
3 points define plane
a*x1+b*y1+c*z1+1=0
a*x2+b*y2+c*z2+1=0
a*x3+b*y2+c*z3+1=0
Solve this system, find a,b,c
Then find z4 for every x4 and y4

When number of base points >3 you can make triangulation for simple and fast 'plain' interpolation or calculate bicubic splines for 'smooth' interpolation
0

LVL 9

Author Comment

ID: 6920249
Hi Mbo,

thanx for comment, but I expect just an implementation code for N points.
Something like Newton - Lagrange polynom solution.

-----
Igor

PS: points can be double.
0

LVL 1

Expert Comment

ID: 6920667
Hi ITugay
Is your point mesh rectangular or irregular?
-----
Boris
0

LVL 1

Expert Comment

ID: 6920671
Hi ITugay
Is your point mesh rectangular or irregular?
-----
Boris
0

LVL 10

Expert Comment

ID: 6922570
Hi ITugay,

Here is a solution. You can make things into nice functions Normalize, DotProduct, CrossProduct but I tested this and it works. I used the gamedev site mentioned in the code. I think the return value should not be integer (and thus not rounded).

Good luck,

Jacco

type
TSFPlane = array[0..2] of TSFPoint;

function GetValueFor(X,Y: Integer; Points: TSFPlane): Integer;
var
dx1, dy1, dv1, dx2, dy2, dv2: Integer;
nx, ny, nv: Integer;
l, rx, ry, rv, d: Extended;
begin
// http://www.gamedev.net/reference/articles/article1443.asp
// check if the three points are not on one line !! (otherwsie this wont work)
dx1 := Points[1].X - Points[0].X;
dy1 := Points[1].Y - Points[0].Y;
dv1 := Points[1].V - Points[0].V;
dx2 := Points[2].X - Points[0].X;
dy2 := Points[2].Y - Points[0].Y;
dv2 := Points[2].V - Points[0].V;
// calcultae cross product (perpendular vector on surface = normal)
nx := dy1 * dv2 - dv1 * dy2;
ny := dv1 * dx2 - dx1 * dv2;
nv := dx1 * dy2 - dy1 * dx2;
// calculate length of normal
l := sqrt(sqr(nx) + sqr(ny) + sqr(nv));
// normalize it (make length = 1)
rx := nx / l;
ry := ny / l;
rv := nv / l;
// calculate distance (intersection of plane with v axis)
d := -(rx * Points[0].X + ry * Points[0].Y + rv * Points[0].V);
// our plane is rx * x + ry * y + rz * z + d = 0 (or z = -rx/rz * x - ry/rz * y - d / rz)
Result := Round( -(rx/rv) * x - (ry/rv) * y - (d/rv) );
end;

0

LVL 9

Author Comment

ID: 6922649
Hi all,

MBo, mesh points irregular.

Jacco, seems it works only for three points, but how about N points?.

Here is an additional explanation. What I'm trying to do.
I have some paper documents and it's images produced by  digital photo camera. The  problem is that images has different lumonicity and color balance in different parts of image. It can be result of flash or wrong light location. I have an utility that makes some preparation of the image for further usage in another application. One of functions should be color balance correction. User can place few points on image where color was white on original. After that, I know location of points and exactly color in that places, it should be possible to correct colors entire image. So, I have to calculate R,G, B correction for all pixels of image.

It's not commercial product, so, I have no time to explore new (for me) branch of science. And hope for some help here. I'm going to give 600 pts with A grade for fully working source code.

-----
Igor.

0

LVL 10

Expert Comment

ID: 6922653
The only thing you need then is a triangulation algorithm.. I' working on it :-).

My idea is:

- triangulate points
- calculate triagle that contains the point
- use previously defined function to calculate value

Regards Jacco
0

LVL 9

Author Comment

ID: 6922672
Hi Jacco,

>> I' working on it :-).
fine :-)

but be carefull, some point can be located outside of any triangle. Let say user clicked only on few points that is not on the edge of image, but area outside of points can be calculated too. Seems you have to calculate edge mesh points before. But for doing this you have to use the same task. Looks like a loop.

----
Igor.

btw. now I'm too trying to build some algorithm, the deal is that I do not need to calculate with high precision, only very approximated values.
0

LVL 9

Author Comment

ID: 6922679
More details,

now I'm trying to make 1D calculations instead of 2D.
step 1. reproject all mesh points to single line relative X coordinte with distance proportionals.
step2. calculate first vlaue for point.
step3. reproject all mesh points to single line relative Y coordinte with distance proportionals.
step 4. calculate second value for point.
step 5. value  = (first + second) div 2.

not sure it will works, just another try :-)

-----
Igor
0

LVL 10

Expert Comment

ID: 6922691
Here is the triangulation module I worked on. I initialized the points list with the four corner points of the image. (Set lighting to 100% on those).

By clicking more points will be added and triangulation is
recalculated.

Is this something we can work with?

unit Unit1;

interface

uses
Windows, Messages, SysUtils, Variants, Classes, Graphics, Controls, Forms,
Dialogs, ExtCtrls, Contnrs, StdCtrls;

type
TVoronoiPoint = class
public
X: Integer;
Y: Integer;
Z: Integer;
constructor Create(aX, aY: Integer);
end;

TVoronoiTriangle = class
p1, p2, p3: TVoronoiPoint;
constructor Create(aP1, aP2, aP3: TVoronoiPoint);
procedure Draw(aCanvas: TCanvas);
end;

TForm1 = class(TForm)
pbxVoronoi: TPaintBox;
Label1: TLabel;
procedure pbxVoronoiPaint(Sender: TObject);
procedure pbxVoronoiMouseDown(Sender: TObject; Button: TMouseButton; Shift: TShiftState; X, Y: Integer);
procedure FormCreate(Sender: TObject);
procedure FormDestroy(Sender: TObject);
private
{ Private declarations }
fPoints: TObjectList;
fTriangles: TObjectList;
public
{ Public declarations }
end;

var
Form1: TForm1;

implementation

{\$R *.dfm}

{ TVoronoiPoint }

constructor TVoronoiPoint.Create(aX, aY: Integer);
begin
X := aX;
Y := aY;
Z := Sqr(X) + Sqr(Y);
end;

{ TVoronoiTriangle }

constructor TVoronoiTriangle.Create(aP1, aP2, aP3: TVoronoiPoint);
begin
p1 := aP1;
p2 := aP2;
p3 := aP3;
end;

procedure TForm1.pbxVoronoiPaint(Sender: TObject);
var
pi, pj, pk, pm: TVoronoiPoint;
xn, yn, zn: Extended;
i, j, k, m, idx: Integer;
begin
if fTriangles.Count = 0 then
begin
// paint the voronoi
for i := 0 to fPoints.Count-1 do
begin
pi := TVoronoiPoint(fPoints[i]);
for j := i + 1 to fPoints.Count-1 do
begin
pj := TVoronoiPoint(fPoints[j]);
for k := i + 1 to fPoints.Count-1 do
begin
pk := TVoronoiPoint(fPoints[k]);
zn := (pj.x - pi.x) * (pk.y - pi.y) - (pk.x - pi.x) * (pj.y - pi.y);
if (j = k) or (zn > 0) then
Continue;
xn := (pj.y - pi.y) * (pk.z - pi.z) - (pk.y - pi.y) * (pj.z - pi.z);
yn := (pk.x - pi.x) * (pj.z - pi.z) - (pj.x - pi.x) * (pk.z - pi.z);
idx := -1;
for m := 0 to fPoints.Count-1 do
begin
pm := TVoronoiPoint(fPoints[m]);
if (m<>i) and (m<>j) and (m<>k) and
(((pm.x - pi.x) * xn + (pm.y - pi.y) * yn + (pm.z - pi.z) * zn) > 0) then
begin
idx := m;
Break;
end;
end;
if idx = -1 then
begin
//for p := 0 to 2 do
//begin

//end;
end;
end;
end;
end;
end;
pbxVoronoi.Canvas.Brush.Color := clWhite;
pbxVoronoi.Canvas.FillRect(pbxVoronoi.ClientRect);
for i := 0 to fTriangles.Count-1 do
begin
TVoronoiTriangle(fTriangles[i]).Draw(pbxVoronoi.Canvas);
end;
end;

procedure TForm1.pbxVoronoiMouseDown(Sender: TObject; Button: TMouseButton; Shift: TShiftState; X, Y: Integer);
begin
fTriangles.Clear;
pbxVoronoi.Invalidate;
end;

procedure TForm1.FormCreate(Sender: TObject);
begin
fPoints := TObjectList.Create;
fTriangles := TObjectList.Create;
end;

procedure TForm1.FormDestroy(Sender: TObject);
begin
fPoints.Free;
fTriangles.Free;
end;

procedure TVoronoiTriangle.Draw(aCanvas: TCanvas);
begin
aCanvas.MoveTo(p1.X, p1.Y);
aCanvas.LineTo(p2.X, p2.Y);
aCanvas.LineTo(p3.X, p3.Y);
aCanvas.LineTo(p1.X, p1.Y);
end;

end.
0

LVL 9

Author Comment

ID: 6922722
Hi Jacco,

seems you on right way. Only one problem, corner lighting  can not be 100%. It should be calculated using at least two points. Let say if lighting tendention from left to right is raised, then right point should be lighter then last right point placed by user. User can't click on some corners because there can be colored zone.

If you'd like to see which images I mean, then let me know, I can send it to you at Monday. Now is slightly later here, so, see you at Monday.

best wishes,
Igor.
0

LVL 10

Expert Comment

ID: 6922764
Ok no edge lighting of 100%

What if I come up with an algorithm that chooses a triangle also if it is outside of all triangles. You can use the GetValueFor for that triangle then.

Let me know what you need.

Regards Jacco
0

LVL 9

Author Comment

ID: 6924545
Hi Jacco,

>> What if I come up with an algorithm that chooses a
>> triangle also if it is outside of all triangles.

yes, you are right, because only three points can be at image.

The way I'm using now:

User makes a picture of empty page on the same place with the same conditions. This image  used to correct rest of images. Stupid way, but works fine :-)
This is why I have to calculate empty page programmatically.

-----
Igor.
0

LVL 9

Author Comment

ID: 6927453
OK,

seems I found a solution:

type

TPoint3D = record
X: Integer;
Y: Integer;
Z: Integer;
end;

TPoints3D = array of TPoint3D;

implementation

function GetZForXY(X, Y: Integer; var PTS: TPoints3D): Integer;
var
D: Real;
I: Integer;
DX, DY: Real;
S1, S2: Real;
begin
S1 := 0;
S2 := 0;
for I := 0 to High(PTS) do
begin
if (X = PTS[I].X) and (Y = PTS[I].Y) then
Result := PTS[I].Z;
Exit;
end;
DX := PTS[I].X - X;
DY := PTS[I].Y - Y;
D := SQRT(DX * DX + DY * DY);

S1 := S1 + PTS[I].Z / D;
S2 := S2 + 1 / D;
end;

Result := Round(S1 / S2);

end;

-------
Igor.
0

LVL 10

Expert Comment

ID: 6929130
This way all lighting points of the image influence every point in the image (with the nearest point having the most influence). Will probably work :)

The solution I am suggesting only uses the three points of the triangle the point lies in (or a triangle on the edge of the mesh if the queried point is outside of all triangles).

Would you like to compare the results?

Regards Jacco
0

LVL 9

Author Comment

ID: 6930408
Hi Jacco,

yes, all points affect to every point. But it is possible to do with only three nearest points, it takes much more code. Result looks the same (at least for color correction). To test this logic I built some gradient bitmaps in PhotoShop. Then place few points on bitmap and make color correction. The result  looking good. The advantage of this way is that possible to place only two or even one point (for linear gradient or image without gradient but wrong "white" color). Also it very simple logic. Seems it cover all my demands.

btw, I'm interesting in another solutions. And don't worry, at least 300 * A points are yours :-)

------
Igor.
0

LVL 10

Accepted Solution

Jacco earned 1200 total points
ID: 6930501
If you solution works it is good :-) It is much simpler. I just made this triangulation and am eager to put it to use somewhere, The triangulation method is known as Delanauy triangulation. I implemeted it to show someone on EE how to implement voronoi diagrams :) (So not much code wasted)

Thanks for the points, although I don't know if I deserve them :)

Regards Jacco
0

LVL 9

Author Comment

ID: 6930536
Hi Jacco,

don't worry, you are really deserved the points. I can imagine a lot of code necessary to implement your way. This is a main reason why I prefer "Inverse Distances Weighting" method.

By the way, if you are going to continue you exploration, please let me know in this thread when it will be finished. Just interesting in it :-)

Regards,

Igor.
0

## Featured Post

Question has a verified solution.

If you are experiencing a similar issue, please ask a related question

In this tutorial I will show you how to use the Windows Speech API in Delphi. I will only cover basic functions such as text to speech and controlling the speed of the speech. SAPI Installation First you need to install the SAPI type library, th…
In my programming career I have only very rarely run into situations where operator overloading would be of any use in my work.  Normally those situations involved math with either overly large numbers (hundreds of thousands of digits or accuracy re…
Visualize your data even better in Access queries. Given a date and a value, this lesson shows how to compare that value with the previous value, calculate the difference, and display a circle if the value is the same, an up triangle if it increased…
In this video, Percona Director of Solution Engineering Jon Tobin discusses the function and features of Percona Server for MongoDB. How Percona can help Percona can help you determine if Percona Server for MongoDB is the right solution for …
###### Suggested Courses
Course of the Month11 days, left to enroll

#### 719 members asked questions and received personalized solutions in the past 7 days.

Join the community of 500,000 technology professionals and ask your questions.