• Status: Solved
  • Priority: Medium
  • Security: Public
  • Views: 383
  • Last Modified:

Interpolation, approximation?

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
task.

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
ITugay
Asked:
ITugay
  • 9
  • 6
  • 3
1 Solution
 
MBoCommented:
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
 
ITugayAuthor Commented:
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
 
MBoCommented:
Hi ITugay
Is your point mesh rectangular or irregular?
-----
Boris
0
Free Tool: SSL Checker

Scans your site and returns information about your SSL implementation and certificate. Helpful for debugging and validating your SSL configuration.

One of a set of tools we are providing to everyone as a way of saying thank you for being a part of the community.

 
MBoCommented:
Hi ITugay
Is your point mesh rectangular or irregular?
-----
Boris
0
 
JaccoCommented:
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
 
ITugayAuthor Commented:
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
 
JaccoCommented:
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
 
ITugayAuthor Commented:
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
 
ITugayAuthor Commented:
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
 
JaccoCommented:
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
            fTriangles.Add(TVoronoiTriangle.Create(pi, pj, pk));
            //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
  // add point
  fTriangles.Clear;
  fPoints.Add(TVoronoiPoint.Create(X,Y));
  pbxVoronoi.Invalidate;
end;

procedure TForm1.FormCreate(Sender: TObject);
begin
  fPoints := TObjectList.Create;
  fTriangles := TObjectList.Create;
  fPoints.Add(TVoronoiPoint.Create(0,0));
  fPoints.Add(TVoronoiPoint.Create(pbxVoronoi.Width-1,0));
  fPoints.Add(TVoronoiPoint.Create(0,pbxVoronoi.Height-1));
  fPoints.Add(TVoronoiPoint.Create(pbxVoronoi.Width-1,pbxVoronoi.Height-1));
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
 
ITugayAuthor Commented:
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
 
JaccoCommented:
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
 
ITugayAuthor Commented:
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
 
ITugayAuthor Commented:
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
 
JaccoCommented:
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
 
ITugayAuthor Commented:
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
 
JaccoCommented:
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
 
ITugayAuthor Commented:
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

2018 Annual Membership Survey

Here at Experts Exchange, we strive to give members the best experience. Help us improve the site by taking this survey today! (Bonus: Be entered to win a great tech prize for participating!)

  • 9
  • 6
  • 3
Tackle projects and never again get stuck behind a technical roadblock.
Join Now