3. Tipps zur Handhabung von Festkommazahlen

Tipps zur Handhabung von Festkommazahlen

  • ?

    Ich würde gerne die Lineare-Algebra-Routinen in meinem Rasterizer mittels Festkommazahlen (int32_t) berechnen.

    Vorherige Versuche Festkommazahlen zu integrieren scheiterten, da

    • Überlauf - der Rasterizer hat dann Dreiecke in Dreiecke gezeichnet (ähnlich wie bei Fraktalen).
    • Ich mir nicht immer sicher bin, wann ich die Skalierung teilen (bei Multiplikation) oder multiplizieren (bei Division) soll.
    • Rundungsfehler (Routine erkennt nicht, dass der Punkt im Dreieck ist und der Pixel gezeichnet werden soll.)

    Habe es zuvor auch mit Wert-Sättigung versucht (also Wert-Grenzen gesetzt), aber zu aufwendig als dass es nützt.

    In Ordnung, nun ans Eingemachte:

    Ich möchte 12 Bits für die Nachkommastellen.
    Somit wäre 2^12 = 4096 die "neue" Nummer 1.
    Dann möchte ich 19 Bits (verschone das 32ste Bit), womit ich dann bei 2^18 - 1 = +262143 und -262144 wäre.
    Wir haben somit eine "Körnung" von 1.0 / 2^12 = 0.000244140625.
    Reicht mir.

    Wenn ich 2 solcher Festkommazahlen multipliziere, dann komme ich mit den Nachkommastellen maximal auf 2^12 * 2^12 = 2^24.
    Bei Division 2^12 / 2^12 = 1.

    Bei der Multiplikation habe ich also schnell eine Überlaufgefahr, oder?

    Soll ich dann auf int64_t konvertieren, da 19 + 24 = 43 und 43 < 64 Bits?
    (Selbst, wenn ich den Vorkommabereich miteinbringe, komme ich auf 19 + 19 + 24 = 62 < 64.)

    static_cast<int32_t>((static_cast<int64_t>(a) * static_cast<int64_t>(b) + 2048) >> 12);

    Ich könnte ja auch die restlichen 7 Bits (32 - 24) verwenden.
    Aber dann wäre mein Vorkommabereich bei +63 und -64.

    Also dann doch der doch recht hässliche Cast da oben?

    Macht es Sinn bei Division zu Runden wie bei der Multiplikation mit 2048?

    Sonst mache ich einfach das

    static_cast<int32_t>((static_cast<int64_t>(a) << 12) / static_cast<int64_t>(b));

    Beide Berechnungen (Multiplikation mit Rundung und Division ohne Rundung) sind jetzt in Funktionen "gepackt"

    • QMultiply(int32_t a, int32_t b);
    • QDivide(int32_t a, int32_t b);

    Soweit so gut.
    Schauen wir uns nun das da unten mal an.

    struct Tuple2 {
    int32_t x;
    int32_t y;
};

constexpr Tuple2 operator-(Tuple2 const& a, Tuple2 const& b) {
    return {a.x - b.x, a.y - b.y};
}

constexpr int Cross(Tuple2 const& a, Tuple2 const& b) {
    return (a.x * b.y - a.y * b.x + kOneHalf) >> kFractionalBits; // kOneHalf sei hier einfach nur 2048 zur Rundung
}

    Bei Cross (2D Kreuzprodukt) war ich mir nicht sicher, ob ich zunächst in Ruhe addieren soll und dann die 12-Bit Rechtsverschiebung oder umgekehrt bzw. einfach 2-mal QMultiply statt *.

    Tuple2 ComputeArealCoordinates(Tuple2 const& p, Tuple2 const& a,
                               Tuple2 const& b, Tuple2 const& c) {
    // Note:
    // p = a + β·(b - a) + γ·(c - a)
    // p - a = β·(b - a) + γ·(c - a)
    // Solve for β, γ using Cramer's rule.
    // Hint: det(a, b) = <e3, a x b> (2D cross product)
    auto const kBMinusA = b - a;
    auto const kCMinusA = c - a;
    auto const kD = Cross(kBMinusA, kCMinusA);

    if (kD < 0) {
        return Tuple2{};
    }

    auto const kPMinusA = p - a;
    auto const kD1 = Cross(kPMinusA, kCMinusA);
    auto const kD2 = Cross(kBMinusA, kPMinusA);
    // auto const kReciprocal = QDivide(kScalingFactor, kD);
    auto const kBeta = QDivide(kD1, kD);
    auto const kGamma = QDivide(kD2, kD);

    return Tuple2{kBeta, kGamma};
}

    Das mit kReciprocal = QDivide(kScalingFactor, kD); und anschließendem QMultiply hat leider nicht geklappt.
    Wollte eine Division sparen ...

    Was meint ihr? Was könnte ich besser machen?

    Ich gehe mal davon aus, dass viele Entwickler hier im Embeddedbereich arbeiten, wo solche Sachen gefragt sind.

    0
  • ?

    Also, ich glaube ich habe die Frage mittlerweile mir selbst beantwortet:

    
#include "rasterizer.h"

#include <limits.h>
#include <stddef.h>

#define ZERO 0
#define ZERO4 \
  { ZERO, ZERO, ZERO, ZERO }

#define COMPUTE_ABSOLUTE_VALUE(x) (((x) < 0) ? -(x) : (x))
#define COMPUTE_DETERMINANT(x0, y0, x1, y1) ((x0) * (y1) - (x1) * (y0))
#define COMPUTE_MINIMUM(x, y) ((x) < (y) ? (x) : (y))
#define COMPUTE_MAXIMUM(x, y) ((x) > (y) ? (x) : (y))

#define FRACTIONAL_BITS 16L
#define SCALED_ONE (1L << FRACTIONAL_BITS)
#define SCALED_ONE_HALF (1L << (FRACTIONAL_BITS - 1L))

#define TO_FIXED_POINT(x) ((long)(x) << FRACTIONAL_BITS)
#define TO_INTEGER(x) ((long)(x) >> FRACTIONAL_BITS)
#define MULTIPLY_AND_SCALE(x, y) (((long)(x) * (long)(y)) >> FRACTIONAL_BITS)
#define DIVIDE_AND_SCALE(x, y) (((long)(x) << FRACTIONAL_BITS) / (long)(y))

/* Rectangle layout: {x_min, x_max, y_min, y_max}. */
static void ComputeBbox(int ax, int ay, int bx, int by, int cx, int cy,
                        int (*rectangle)[4]) {
  int* ptr = *rectangle;
  /* Compute x min/max and y min/max. */
  ptr[0] = COMPUTE_MINIMUM(ax, COMPUTE_MINIMUM(bx, cx));
  ptr[1] = COMPUTE_MAXIMUM(ax, COMPUTE_MAXIMUM(bx, cx));
  ptr[2] = COMPUTE_MINIMUM(ay, COMPUTE_MINIMUM(by, cy));
  ptr[3] = COMPUTE_MAXIMUM(ay, COMPUTE_MAXIMUM(by, cy));
}

void RasterizeTriangle(int ax, int ay, int bx, int by, int cx, int cy,
                       ComputeColorFunction* compute_color,
                       DrawPixelWithColorFunction* draw_pixel, void* data) {
  /* W1 := OB - OA */
  int const kW1X = bx - ax;
  int const kW1Y = by - ay;

  /* W2 := OC - OA */
  int const kW2X = cx - ax;
  int const kW2Y = cy - ay;

  /* Det(w1, w2) = <w1 x w2, e3> = |w1| * |w2| * sin(α) */
  int const kD = COMPUTE_DETERMINANT(kW1X, kW1Y, kW2X, kW2Y);

  /* Cue: |w1| * |w2| * sin(α) */
  int const kDegenerate = kD == 0;

  int rectangle[4] = ZERO4;

  /* Ignore degenerate triangles. */
  if (kDegenerate == 1) {
    return;
  }

  ComputeBbox(ax, ay, bx, by, cx, cy, &rectangle);

  {
    /* Minimum bounding rectangle. */
    int const kMinX = rectangle[0];
    int const kMaxX = rectangle[1];
    int const kMinY = rectangle[2];
    int const kMaxY = rectangle[3];

    /* Cue: float kDReciprocal = 1.0f / (float)kD; */
    long const kReciprocal = DIVIDE_AND_SCALE(SCALED_ONE, TO_FIXED_POINT(kD));

    int y = 0;
    int x = 0;
    for (y = kMinY; y < kMaxY; ++y) {
      for (x = kMinX; x < kMaxX; ++x) {
        /* P - V0 */
        int const kWX = x - ax;
        int const kWY = y - ay;

        /* Det(w, w2) = <w x w2, e3> = |w| * |w2| * sin(α) */
        int const kD1 = COMPUTE_DETERMINANT(kWX, kWY, kW2X, kW2Y);

        /* Det(w1, w) = <w1 x w, e3> = |w1| * |w| * sin(α) */
        int const kD2 = COMPUTE_DETERMINANT(kW1X, kW1Y, kWX, kWY);

        /* kS1 := D1 / D */
        long const kS1 = MULTIPLY_AND_SCALE(TO_FIXED_POINT(kD1), kReciprocal);

        /* kS2 := D2 / D */
        long const kS2 = MULTIPLY_AND_SCALE(TO_FIXED_POINT(kD2), kReciprocal);

        /* kS3 := 1 - S1 - S2 */
        long const kS3 = SCALED_ONE - kS1 - kS2;

        /* Hint: (kS1 | kS2 | kS3) >= 0 ~ kS1 >= 0 && kS2 >= 0 && kS3 >= 0 */
        int const kPointInside = (kS1 | kS2 | kS3) >= 0L;

        if (kPointInside == 0) {
          continue;
        }

        /* Draw pixel. */
        {
          unsigned char colors[3] = {0, 0, 0};
          unsigned char* color = colors;
          compute_color(kS1, kS2, kS3, FRACTIONAL_BITS, &colors, data);
          draw_pixel(x, y, color[0], color[1], color[2]);
        }
      }
    }
  }
}

    Mit diesem Code von mir, kann man, so weit ich in Erfahrung bringen konnte, das 3D-Objkekt dasrstellen.
    Nur hakt es beim Interpolieren mit Texturen.
    Ich glaube, der Fehler ist beim Interpolieren mit den Festkommawerten von kS1 (alpha), kS2 (beta) und kS3 (gamma) zugreife und dann mit FRACTIONAL_BITS runterskaliere auf float.
    Nebenbei angemerkt, kS1, kS2 kS3 sind schreckliche Namen.
    Wäre vllt. besser, diese so zu nennen:

    • kScaledAlpha
    • kScaledBeta
    • kScaledGamma

    Edit: Und das Problem war, dass ich zuvor 12-Bits verwendet habe statt mindestens 16-Bit.
    Genauigkeit ist wichtig.
    Das war der Hauptgrund denke ich.

    0
Anmelden zum Antworten