Autor Wątek: Kolizja AABB i promienia  (Przeczytany 1602 razy)

Offline Witek9002

  • Użytkownik

# Sierpień 30, 2009, 15:36:52
Zna ktoś może dobry algorytm znajdywania kolizjii między promieniem i AABB w 3D, który zwraca również punkt, gdzie nastąpiła kolizja? Próbowałem tego: http://jgt.akpeters.com/papers/MahovskyWyvill04/JGT-float.zip Działa bardzo szybko, tylko że punkt kolizji jest źle obliczany...     :-\

Aktualnie muszę uzywać tego:
function RayAABBIntersection(const Ray: PRay; const AABB: PAABB; var t: Single) : Boolean;
var
  tnear, tfar, t1, t2, temp: Single;
begin
    tnear := -1e6;
    tfar := 1e6;

    t1 := (AABB^.Min[0] - Ray^.Origin[0]) / Ray^.Dir[0];
    t2 := (AABB^.Max[0] - Ray^.Origin[0]) / Ray^.Dir[0];

    if (t1 > t2) then
    begin
        temp := t1;
        t1 := t2;
        t2 := temp;
    end;

    if (t1 > tnear) then
        tnear := t1;
    if (t2 < tfar)then
        tfar := t2;

    if (tnear > tfar) then
    begin
        Result := false;
        exit;
    end;

    if (tfar < 0.0) then
    begin
        Result := false;
        exit;
    end;


    t1 := (AABB^.Min[1] - Ray^.Origin[1]) / Ray^.Dir[1];
    t2 := (AABB^.Max[1] - Ray^.Origin[1]) / Ray^.Dir[1];

    if (t1 > t2) then
    begin
        temp := t1;
        t1 := t2;
        t2 := temp;
    end;

    if (t1 > tnear) then
        tnear := t1;
    if (t2 < tfar)then
        tfar := t2;

    if (tnear > tfar) then
    begin
        Result := false;
        exit;
    end;

    if (tfar < 0.0) then
    begin
        Result := false;
        exit;
    end;

    t1 := (AABB^.Min[2] - Ray^.Origin[2]) / Ray^.Dir[2];
    t2 := (AABB^.Max[2] - Ray^.Origin[2]) / Ray^.Dir[2];

    if (t1 > t2) then
    begin
        temp := t1;
        t1 := t2;
        t2 := temp;
    end;

    if (t1 > tnear) then
        tnear := t1;
    if (t2 < tfar)then
        tfar := t2;

    if (tnear > tfar) then
    begin
        Result := false;
        exit;
    end;

    if (tfar < 0.0) then
    begin
        Result := false;
        exit;
    end;

    t := tnear;
    Result := true;
end;

Działa w 100%, ale jest strasznie wolne (okolo 4x od tych w JGT-float.zip)

Z góry dziękuję za pomoc
« Ostatnia zmiana: Sierpień 30, 2009, 15:45:57 wysłana przez Witek9002 »

Offline Mr. Spam

  • Miłośnik przetworów mięsnych

Offline yarpen

  • Użytkownik

# Sierpień 30, 2009, 16:00:38
Google: Kay-Kajiya slab volume

Offline Reg

  • Administrator
    • Adam Sawicki - Home Page

# Sierpień 30, 2009, 21:18:16
Ja używam takiej funkcji. Już nie pamiętam skąd ją mam :)

Dostając parametr t (który jest wyrażony w wielokrotnościach długości RayDir, bo ten nie musi być znormalizowany) można wyliczyć punkt przecięcia ze wzoru:
Punkt_przeciecia = RayOrig + RayDir*t;

bool RayToBox(float *OutT, const VEC3 &RayOrig, const VEC3 &RayDir, const BOX &Box)
{
bool inside = true;
float xt, xn;

if (RayOrig.x < Box.Min.x)
{
xt = Box.Min.x - RayOrig.x;
xt /= RayDir.x;
xn = -1.0f;
inside = false;
}
else if (RayOrig.x > Box.Max.x)
{
xt = Box.Max.x - RayOrig.x;
xt /= RayDir.x;
xn = 1.0f;
inside = false;
}
else
xt = -1.0f;

float yt, yn;

if (RayOrig.y < Box.Min.y)
{
yt = Box.Min.y - RayOrig.y;
yt /= RayDir.y;
yn = -1.0f;
inside = false;
}
else if (RayOrig.y > Box.Max.y)
{
yt = Box.Max.y - RayOrig.y;
yt /= RayDir.y;
yn = 1.0f;
inside = false;
}
else
yt = -1.0f;

float zt, zn;

if (RayOrig.z < Box.Min.z)
{
zt = Box.Min.z - RayOrig.z;
zt /= RayDir.z;
zn = -1.0f;
inside = false;
}
else if (RayOrig.z > Box.Max.z)
{
zt = Box.Max.z - RayOrig.z;
zt /= RayDir.z;
zn = 1.0f;
inside = false;
}
else
zt = -1.0f;

if (inside)
{
*OutT = 0.0f;
return true;
}

// Select the farthest plane - this is the plane of intersection
int plane = 0;

float t = xt;
if (yt > t)
{
plane = 1;
t = yt;
}

if (zt > t)
{
plane = 2;
t = zt;
}

// Check if the point of intersection lays within the box face

switch(plane)
{
case 0: // ray intersects with yz plane
{
float y = RayOrig.y + RayDir.y * t;
if (y < Box.Min.y || y > Box.Max.y) return false;
float z = RayOrig.z + RayDir.z * t;
if (z < Box.Min.z || z > Box.Max.z) return false;
}
break;
case 1: // ray intersects with xz plane
{
float x = RayOrig.x + RayDir.x * t;
if (x < Box.Min.x || x > Box.Max.x) return false;
float z = RayOrig.z + RayDir.z * t;
if (z < Box.Min.z || z > Box.Max.z) return false;
}
break;
default:
case 2: // ray intersects with xy plane
{
float x = RayOrig.x + RayDir.x * t;
if (x < Box.Min.x || x > Box.Max.x) return false;
float y = RayOrig.y + RayDir.y * t;
if (y < Box.Min.y || y > Box.Max.y) return false;
}
break;
}

*OutT = t;
return true;
}

Offline Witek9002

  • Użytkownik

# Sierpień 30, 2009, 21:33:17
Cytuj
Dostając parametr t (który jest wyrażony w wielokrotnościach długości RayDir, bo ten nie musi być znormalizowany) można wyliczyć punkt przecięcia ze wzoru:
Punkt_przeciecia = RayOrig + RayDir*t;

Wiem o tym. Poza tym już sobie poradziłem, znalazłem błąd, wszystko już działa OK. Algorytm jest strasznie dlugi, ale bardzo szybki, raczej szybszy od twojego, Reg. Jak ktoś chce, to mogę go zamieścić

Offline Reg

  • Administrator
    • Adam Sawicki - Home Page

# Sierpień 31, 2009, 21:25:07
Ja chcę! Zamieść :)

Offline Witek9002

  • Użytkownik

# Sierpień 31, 2009, 22:28:19
Bardzo proszę.  ;)
Niestety jest w Delphi, ale mam nadzieję, że sobie poradzisz.
PS. Sorry za bałagan z akapitami.

Przed wywołaniem RayAABBIntersection(@Ray, @Box, t) trzeba użyć RayClassificate(@Ray);

Kod:

// deklaracja typow itd...
type
  TVector3f = array[0..2] of single;

  PRay = ^TRay;
  TRay = record
    Origin, Dir: TVector3f;
    classification: Byte;
  end;

  PAABB = ^TAABB;
  TAABB =
    record
      Min, Max: TVector3f;
    end;


const
  RAY_MMM = 0;
  RAY_MMP = 1;
  RAY_MPM = 2;
  RAY_MPP = 3;
  RAY_PMM = 4;
  RAY_PMP = 5;
  RAY_PPM = 6;
  RAY_PPP = 7;




//wykrywanie kolizji AABB-promien
function RayAABBIntersection(const Ray: PRay; const B: PAABB; var t: Single) : Boolean;
var
  xa, xb, ya, yb, za, zb, t1, t2: Single;
begin
  case Ray^.classification of
    RAY_MMM:
      begin
if ((Ray^.Origin[0] < B^.Min[0]) or (Ray^.Origin[1] < B^.Min[1]) or (Ray^.Origin[2] < B^.Min[2])) then
        begin
  Result := false;
            exit;
        end;

xa := B^.Min[0] - Ray^.Origin[0];
ya := B^.Min[1] - Ray^.Origin[1];
za := B^.Min[2] - Ray^.Origin[2];
xb := B^.Max[0] - Ray^.Origin[0];
yb := B^.Max[1] - Ray^.Origin[1];
zb := B^.Max[2] - Ray^.Origin[2];

if (Ray^.Dir[0] * ya - Ray^.Dir[1] * xb < 0) or
(Ray^.Dir[0] * yb - Ray^.Dir[1] * xa > 0) or
(Ray^.Dir[0] * zb - Ray^.Dir[2] * xa > 0) or
(Ray^.Dir[0] * za - Ray^.Dir[2] * xb < 0) or
(Ray^.Dir[1] * za - Ray^.Dir[2] * yb < 0) or
(Ray^.Dir[1] * zb - Ray^.Dir[2] * ya > 0) then
          Result := false
        else
        begin
            t := xb / Ray^.Dir[0];
      t1 := yb / Ray^.Dir[1];
      if (t1 > t) then
      t := t1;
      t2 := zb / Ray^.Dir[2];
      if (t2 > t) then
      t := t2;

            Result := true;
        end;
      end;

    RAY_MMP:
      begin

// case MMP: side(R,HD) < 0 or side(R,FB) > 0 or side(R,HG) > 0 or side(R,AB) < 0 or side(R,DA) < 0 or side(R,GF) > 0 to miss

if ((Ray^.Origin[0] < B^.Min[0]) or (Ray^.Origin[1] < B^.Min[1]) or (Ray^.Origin[2] > B^.Max[2])) then
        begin
  Result := false;
            exit;
        end;

xa := B^.Min[0] - Ray^.Origin[0];
ya := B^.Min[1] - Ray^.Origin[1];
za := B^.Min[2] - Ray^.Origin[2];
xb := B^.Max[0] - Ray^.Origin[0];
yb := B^.Max[1] - Ray^.Origin[1];
zb := B^.Max[2] - Ray^.Origin[2];

if (Ray^.Dir[0] * ya - Ray^.Dir[1] * xb < 0) or
(Ray^.Dir[0] * yb - Ray^.Dir[1] * xa > 0) or
(Ray^.Dir[0] * zb - Ray^.Dir[2] * xb > 0) or
(Ray^.Dir[0] * za - Ray^.Dir[2] * xa < 0) or
(Ray^.Dir[1] * za - Ray^.Dir[2] * ya < 0) or
(Ray^.Dir[1] * zb - Ray^.Dir[2] * yb > 0) then
          Result := false
        else
        begin
            t := xb / Ray^.Dir[0];
      t1 := yb / Ray^.Dir[1];
      if (t1 > t) then
      t := t1;
      t2 := za / Ray^.Dir[2];
      if (t2 > t) then
      t := t2;

            Result := true;
        end;
      end;

    RAY_MPM:
      begin
// case MPM: side(R,EA) < 0 or side(R,GC) > 0 or side(R,EF) > 0 or side(R,DC) < 0 or side(R,GF) < 0 or side(R,DA) > 0 to miss

if ((Ray^.Origin[0] < B^.Min[0]) or (Ray^.Origin[1] > B^.Max[1]) or (Ray^.Origin[2] < B^.Min[2])) then
        begin
  Result := false;
            exit;
        end;

xa := B^.Min[0] - Ray^.Origin[0];
ya := B^.Min[1] - Ray^.Origin[1];
za := B^.Min[2] - Ray^.Origin[2];
xb := B^.Max[0] - Ray^.Origin[0];
yb := B^.Max[1] - Ray^.Origin[1];
zb := B^.Max[2] - Ray^.Origin[2];

if (Ray^.Dir[0] * ya - Ray^.Dir[1] * xa < 0) or
(Ray^.Dir[0] * yb - Ray^.Dir[1] * xb > 0) or
(Ray^.Dir[0] * zb - Ray^.Dir[2] * xa > 0) or
(Ray^.Dir[0] * za - Ray^.Dir[2] * xb < 0) or
(Ray^.Dir[1] * zb - Ray^.Dir[2] * yb < 0) or
(Ray^.Dir[1] * za - Ray^.Dir[2] * ya > 0) then
          Result := false
        else
        begin
            t := xb / Ray^.Dir[0];
      t1 := ya / Ray^.Dir[1];
      if (t1 > t) then
      t := t1;
      t2 := zb / Ray^.Dir[2];
      if (t2 > t) then
      t := t2;

            Result := true;
        end;
      end;

    RAY_MPP:
      begin
// case MPP: side(R,EA) < 0 or side(R,GC) > 0 or side(R,HG) > 0 or side(R,AB) < 0 or side(R,HE) < 0 or side(R,CB) > 0 to miss

if ((Ray^.Origin[0] < B^.Min[0]) or (Ray^.Origin[1] > B^.Max[1]) or (Ray^.Origin[2] > B^.Max[2])) then
        begin
  Result := false;
            exit;
        end;

xa := B^.Min[0] - Ray^.Origin[0];
ya := B^.Min[1] - Ray^.Origin[1];
za := B^.Min[2] - Ray^.Origin[2];
xb := B^.Max[0] - Ray^.Origin[0];
yb := B^.Max[1] - Ray^.Origin[1];
zb := B^.Max[2] - Ray^.Origin[2];

if (Ray^.Dir[0] * ya - Ray^.Dir[1] * xa < 0) or
(Ray^.Dir[0] * yb - Ray^.Dir[1] * xb > 0) or
(Ray^.Dir[0] * zb - Ray^.Dir[2] * xb > 0) or
(Ray^.Dir[0] * za - Ray^.Dir[2] * xa < 0) or
(Ray^.Dir[1] * zb - Ray^.Dir[2] * ya < 0) or
(Ray^.Dir[1] * za - Ray^.Dir[2] * yb > 0) then
          Result := false
        else
        begin
            t := xb / Ray^.Dir[0];
      t1 := ya / Ray^.Dir[1];
      if (t1 > t) then
      t := t1;
      t2 := za / Ray^.Dir[2];
      if (t2 > t) then
      t := t2;

            Result := true;
        end;
      end;

    RAY_PMM:
      begin
// case PMM: side(R,GC) < 0 or side(R,EA) > 0 or side(R,AB) > 0 or side(R,HG) < 0 or side(R,CB) < 0 or side(R,HE) > 0 to miss

if ((Ray^.Origin[0] > B^.Max[0]) or (Ray^.Origin[1] < B^.Min[1]) or (Ray^.Origin[2] < B^.Min[2]))  then
        begin
  Result := false;
            exit;
        end;

xa := B^.Min[0] - Ray^.Origin[0];
ya := B^.Min[1] - Ray^.Origin[1];
za := B^.Min[2] - Ray^.Origin[2];
xb := B^.Max[0] - Ray^.Origin[0];
yb := B^.Max[1] - Ray^.Origin[1];
zb := B^.Max[2] - Ray^.Origin[2];

if (Ray^.Dir[0] * yb - Ray^.Dir[1] * xb < 0) or
(Ray^.Dir[0] * ya - Ray^.Dir[1] * xa > 0) or
(Ray^.Dir[0] * za - Ray^.Dir[2] * xa > 0) or
(Ray^.Dir[0] * zb - Ray^.Dir[2] * xb < 0) or
(Ray^.Dir[1] * za - Ray^.Dir[2] * yb < 0) or
(Ray^.Dir[1] * zb - Ray^.Dir[2] * ya > 0) then
          Result := false
        else
        begin
            t := xa / Ray^.Dir[0];
      t1 := yb / Ray^.Dir[1];
      if (t1 > t) then
      t := t1;
      t2 := zb / Ray^.Dir[2];
      if (t2 > t) then
      t := t2;

            Result := true;
        end;
      end;

    RAY_PMP:
      begin
// case PMP: side(R,GC) < 0 or side(R,EA) > 0 or side(R,DC) > 0 or side(R,EF) < 0 or side(R,DA) < 0 or side(R,GF) > 0 to miss

if ((Ray^.Origin[0] > B^.Max[0]) or (Ray^.Origin[1] < B^.Min[1]) or (Ray^.Origin[2] > B^.Max[2])) then
        begin
  Result := false;
            exit;
        end;

xa := B^.Min[0] - Ray^.Origin[0];
ya := B^.Min[1] - Ray^.Origin[1];
za := B^.Min[2] - Ray^.Origin[2];
xb := B^.Max[0] - Ray^.Origin[0];
yb := B^.Max[1] - Ray^.Origin[1];
zb := B^.Max[2] - Ray^.Origin[2];

if (Ray^.Dir[0] * yb - Ray^.Dir[1] * xb < 0) or
(Ray^.Dir[0] * ya - Ray^.Dir[1] * xa > 0) or
(Ray^.Dir[0] * za - Ray^.Dir[2] * xb > 0) or
(Ray^.Dir[0] * zb - Ray^.Dir[2] * xa < 0) or
(Ray^.Dir[1] * za - Ray^.Dir[2] * ya < 0) or
(Ray^.Dir[1] * zb - Ray^.Dir[2] * yb > 0) then
          Result := false
        else
        begin
            t := xa / Ray^.Dir[0];
      t1 := yb / Ray^.Dir[1];
      if (t1 > t) then
      t := t1;
      t2 := za / Ray^.Dir[2];
      if (t2 > t) then
      t := t2;

            Result := true;
        end;

      end;

    RAY_PPM:
      begin
// case PPM: side(R,FB) < 0 or side(R,HD) > 0 or side(R,AB) > 0 or side(R,HG) < 0 or side(R,GF) < 0 or side(R,DA) > 0 to miss

if ((Ray^.Origin[0] > B^.Max[0]) or (Ray^.Origin[1] > B^.Max[1]) or (Ray^.Origin[2] < B^.Min[2])) then
        begin
  Result := false;
            exit;
        end;

xa := B^.Min[0] - Ray^.Origin[0];
ya := B^.Min[1] - Ray^.Origin[1];
za := B^.Min[2] - Ray^.Origin[2];
xb := B^.Max[0] - Ray^.Origin[0];
yb := B^.Max[1] - Ray^.Origin[1];
zb := B^.Max[2] - Ray^.Origin[2];

if (Ray^.Dir[0] * yb - Ray^.Dir[1] * xa < 0) or
(Ray^.Dir[0] * ya - Ray^.Dir[1] * xb > 0) or
(Ray^.Dir[0] * za - Ray^.Dir[2] * xa > 0) or
(Ray^.Dir[0] * zb - Ray^.Dir[2] * xb < 0) or
(Ray^.Dir[1] * zb - Ray^.Dir[2] * yb < 0) or
(Ray^.Dir[1] * za - Ray^.Dir[2] * ya > 0) then
          Result := false
        else
        begin
            t := xa / Ray^.Dir[0];
      t1 := ya / Ray^.Dir[1];
      if (t1 > t) then
      t := t1;
      t2 := zb / Ray^.Dir[2];
      if (t2 > t) then
      t := t2;

            Result := true;
        end;
      end;

    RAY_PPP:
      begin
// case PPP: side(R,FB) < 0 or side(R,HD) > 0 or side(R,DC) > 0 or side(R,EF) < 0 or side(R,HE) < 0 or side(R,CB) > 0 to miss

if ((Ray^.Origin[0] > B^.Max[0]) or (Ray^.Origin[1] > B^.Max[1]) or (Ray^.Origin[2] > B^.Max[2]))  then
        begin
  Result := false;
            exit;
        end;

xa := B^.Min[0] - Ray^.Origin[0];
ya := B^.Min[1] - Ray^.Origin[1];
za := B^.Min[2] - Ray^.Origin[2];
xb := B^.Max[0] - Ray^.Origin[0];
yb := B^.Max[1] - Ray^.Origin[1];
zb := B^.Max[2] - Ray^.Origin[2];

if (Ray^.Dir[0] * yb - Ray^.Dir[1] * xa < 0) or
(Ray^.Dir[0] * ya - Ray^.Dir[1] * xb > 0) or
(Ray^.Dir[0] * za - Ray^.Dir[2] * xb > 0) or
(Ray^.Dir[0] * zb - Ray^.Dir[2] * xa < 0) or
(Ray^.Dir[1] * zb - Ray^.Dir[2] * ya < 0) or
(Ray^.Dir[1] * za - Ray^.Dir[2] * yb > 0) then
          Result := false
        else
        begin
            t := xa / Ray^.Dir[0];
      t1 := ya / Ray^.Dir[1];
      if (t1 > t) then
      t := t1;
      t2 := za / Ray^.Dir[2];
      if (t2 > t) then
      t := t2;

            Result := true;
        end;
      end;
    else
      Result := false;
    end;
end;





// "klasyfikowanie" promienia w zależności od kierunku
procedure RayClassificate(const Ray: PRay);
begin
if(Ray^.Dir[0] < 0) then
begin
if(Ray^.Dir[1] < 0) then
begin
if(Ray^.Dir[2] < 0) then
Ray^.classification := RAY_MMM
else
Ray^.classification := RAY_MMP;
end
else
begin
if(Ray^.Dir[2] < 0) then
Ray^.classification := RAY_MPM
else
Ray^.classification := RAY_MPP;
end
end
else
begin
if(Ray^.Dir[1] < 0) then
begin
if(Ray^.Dir[2] < 0) then
Ray^.classification := RAY_PMM
else
Ray^.classification := RAY_PMP;
end
else
begin
if(Ray^.Dir[2] < 0) then
Ray^.classification := RAY_PPM
else
Ray^.classification := RAY_PPP;
end
end
end;


Offline Witek9002

  • Użytkownik

# Czerwiec 02, 2011, 00:58:48
Tak, wiem że odgrzewam kotlet, ale chciałbym się pochwalić moją nową funkcją do testowania kolizji promienia i AABB. Opiera się na SSE, jest ultra-szybka i zajmuje 27 instrukcji procesora bez żadnego branchingu.
Postanowiłem, że ją tutaj zamieszczę - może przyda się przyszłym pokoleniom ;)

Oto ona:
//struktura opisująca AABB
struct Box
{
__m128 Min;
__m128 Max;
};

//struktura opisująca promień
struct Ray
{
__m128 Origin; //początek
__m128 Dir;  //kierunek
__m128 invDir; //odwrotność Dir
};

inline int IntersectRayBox(const Box &box, const Ray &ray, __m128* pNear, __m128* pFar)
{
__m128 tmp1 = _mm_mul_ps(_mm_sub_ps(box.Min, ray.Origin), ray.invDir);
__m128 tmp2 = _mm_mul_ps(_mm_sub_ps(box.Max, ray.Origin), ray.invDir);
__m128 lmin = _mm_min_ps(tmp1, tmp2);
__m128 lmax = _mm_max_ps(tmp1, tmp2);

__m128 lx = _mm_shuffle_ps(lmin, lmax, 0x0);
__m128 ly = _mm_shuffle_ps(lmin, lmax, 0x55);
__m128 lz = _mm_shuffle_ps(lmin, lmax, 0xAA);

lmin = _mm_max_ps(lx, _mm_max_ps(ly, lz));
lmax = _mm_min_ps(lx, _mm_min_ps(ly, lz));
lmax = _mm_shuffle_ps(lmax, lmax, 0xAA);
lmin = _mm_shuffle_ps(lmin, lmin, 0x55);
*pNear = lmin;
*pFar = lmax;

tmp1 = _mm_cmpge_ps(lmax, _mm_setzero_ps());
tmp2 = _mm_cmpge_ps(lmax, lmin);
return (_mm_movemask_ps(_mm_and_ps(tmp1, tmp2)) & 7) == 7;
}

Funkcja zwraca 0 lub 1. Działa również dla przypadku, w którym początek promienia znajduje się wewnątrz AABB. Poprzez parametry pNear i pFar można pobrać odległość do kolizji z przednią i tylną ścianką.

Uwaga: Należy podać odwrotność znormalizowanego kierunku promienia w ray.invDir!

PS. gdy zobaczyłem funkcje z poprzednich postów, to się roześmiałem :)

Offline Krzysiek K.

  • Redaktor
    • DevKK.net

# Czerwiec 02, 2011, 11:51:27
Cytuj
Tak, wiem że odgrzewam kotlet, ale chciałbym się pochwalić moją nową funkcją do testowania kolizji promienia i AABB. Opiera się na SSE, jest ultra-szybka i zajmuje 27 instrukcji procesora bez żadnego branchingu.
Postanowiłem, że ją tutaj zamieszczę - może przyda się przyszłym pokoleniom ;)
Używałem właściwie dokładnie tego samego algorytmu, tyle że odpowiednik lmin/lmax wyznaczałem używając czegoś w stylu "invdir>0?bmin:bmax", więc min/max na CPU w tym miejscu będzie faktycznie szybsze (no i poza tym Twoja implementacja jest SSE). Ani jedno ani drugie nie zmienia jednak za bardzo wydajności z mojej strony, bo obecnie wykorzystuję implementacje głównie na GPU i w Javie.