Solved

Interpolation, approximation?

Posted on 2002-04-04
18
368 Views
Last Modified: 2010-04-04
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
Comment
Question by:ITugay
  • 9
  • 6
  • 3
18 Comments
 
LVL 1

Expert Comment

by:MBo
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

by:ITugay
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

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

Expert Comment

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

Expert Comment

by:Jacco
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

by:ITugay
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

by:Jacco
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

by:ITugay
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

by:ITugay
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
How your wiki can always stay up-to-date

Quip doubles as a “living” wiki and a project management tool that evolves with your organization. As you finish projects in Quip, the work remains, easily accessible to all team members, new and old.
- Increase transparency
- Onboard new hires faster
- Access from mobile/offline

 
LVL 10

Expert Comment

by:Jacco
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
            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
 
LVL 9

Author Comment

by:ITugay
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

by:Jacco
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

by:ITugay
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

by:ITugay
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

by:Jacco
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

by:ITugay
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

by:
Jacco earned 300 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

by:ITugay
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

What Should I Do With This Threat Intelligence?

Are you wondering if you actually need threat intelligence? The answer is yes. We explain the basics for creating useful threat intelligence.

Join & Write a Comment

Suggested Solutions

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…
Hello everybody This Article will show you how to validate number with TEdit control, What's the TEdit control? TEdit is a standard Windows edit control on a form, it allows to user to write, read and copy/paste single line of text. Usua…
Sending a Secure fax is easy with eFax Corporate (http://www.enterprise.efax.com). First, Just open a new email message.  In the To field, type your recipient's fax number @efaxsend.com. You can even send a secure international fax — just include t…
This tutorial demonstrates a quick way of adding group price to multiple Magento products.

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

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

Join & Ask a Question

Need Help in Real-Time?

Connect with top rated Experts

12 Experts available now in Live!

Get 1:1 Help Now