Υπολογισμός όγκου Simplex σε C++

[H_ANARXIA_EINAI_PSEMA]

Το παρακάτω πρόγραμμα παίρνει ως είσοδο μια συλλογή από σημεία (5-άδες στο συγκεκριμένο πρόγραμμα, αλλάξτε το d) και βρίσκει τον όγκο του simplex που φτιάχνουν τα σημεία αυτά. Αν τυχόν η συλλογή δεν είναι καλή (δηλ. δεν φτιάχνει simplex), το πρόγραμμα σας το λέει.

 

 

ο κώδικας:

 

#include <iostream>
#include <vector>
#include <array>
#include <algorithm>
#include <functional>
#include <cmath>

/* This is the dimension R^d */
constexpr int d = 5;

using namespace std;

/* Test whether the vector v is all zero */
bool zero(const array<double, d>& v)
{
  for(auto x : v) if(x) return false;
  return true;
}

/* Test whether the set of vectors
   is linearly dependent */
bool lindep(vector<array<double, d>> v)
{
  if(v.size() == 0) return false;
  if(v.size() > d) return true;
  for(auto it = v.begin(); it != v.end(); ++it) {
      if(zero(*it)) return true;
      // Find leading coefficient
      int k =
        find_if(it->begin(), it->end(),
                [](double f)->bool{ return f!=0; })
        - it->begin();
      // Zeroes below leading coefficients
      for(auto tmp = it + 1; tmp != v.end(); ++tmp) {
        double c = tmp[0][k]/it[0][k];
        for(int i = k; i < it->size(); i++)
          tmp[0][i] -= c*it[0][i];
      }
    }
  return false;
}

/* Calculate the determinant of a d x d matrix
   Uses the Leibniz formula */
double det(const vector<array<double, d>>& m)
{
  array<double, d> i;
  double rv = 0.0;
  bool sigma = 1;
  for(int j = 0; j < i.size(); j++)
    i[j] = j;
  do {
    double tmp = 1.0;
    for(int j = 0; j < i.size(); j++)
      tmp *= m[j][i[j]];
    rv += (1-2*sigma)*tmp;
    sigma = !sigma;
  } while(next_permutation(i.begin(), i.end()));
  return rv;
}

int main(void)
{
  vector<array<double, d>> v, w; /* collections of points */
  array<double, d> tmp;
  string s = "";
  cout << "Enter points in format n1 ... n"
       << d << " separated by newline:" << endl
       << "(give " << d << "+1 points to "
       << "test for simplex)" << endl;
  /* Read points until end of file */
  while(1) {
    for(auto& x : tmp) cin >> x;
    if(cin.eof()) break;
    v.push_back(tmp);
  }
  if(v.empty()) s = "in";
  else {
    w = v;
    tmp = w.back();
    w.pop_back();
    if(!(zero(tmp) && w.empty())) {
      for(auto& x : w)
        for(int i = 0; i < x.size(); i++)
          x[i] -= tmp[i];
      if(!lindep(w)) s = "in";
    }
  }
  if(v.size() != d + 1)
    cout << "Affinely " << s << "dependent." << endl;
  else if(s == "in") {
    function<double(double)> f = [&f](double n){return n?n*f(n-1):1;};
    cout << "Simplex volume: " << abs(det(w))/f(d) << endl;
  } else cout << "Not a simplex." << endl;
  return 0;
}

 

Παράδειγμα

 

anarchy@bitch % ./a.out
Enter points in format n1 ... n5 separated by newline:
(give 5+1 points to test for simplex)
0 0 0 0 0
1 0 0 0 0
0 0 1 0 0
0 1 0 0 0
0 0 0 1 0
0 0 0 0 1
Simplex volume: 0.00833333

[V.I.Smirnov]

Ενδιαφέρον, αν και δεν εντυπωσιάζομαι ως προς την ιδέα.

Η επιφάνεια και ο όγκος οποιουδήποτε πολυγώνου ή πολυέδρου μπορούν να βρεθούν πολύ εύκολα με

τα θ. Stokes και Gauss.

Είναι ένα άθροισμα εξωτερικών γινομένων των διανυσμάτων θέσης των κορυφών και των καθέτων των πλευρών.

Θυμάμαι μάλιστα ένα σύντομο άρθρο επ' αυτού σε κάποιον τόμο των graphics gems (1991).

 

-

[παπι]

Όλο λέω να γράψω ένα πρόγραμμα που προβάλει 4d και όλο δεν το γράφω...

 

N1 για αυτό!που έγραψες αν και δεν ξέρω τι είναι.

[H_ANARXIA_EINAI_PSEMA]

Όλο λέω να γράψω ένα πρόγραμμα που προβάλει 4d και όλο δεν το γράφω...

 

N1 για αυτό!που έγραψες αν και δεν ξέρω τι είναι.

Παίρνωντας διανύσματα (σημεία) στο χώρο, μπορείς να σχηματίσεις την κλειστή θήκη δηλαδή αυτό που θα σχηματιζόταν αν προσπαθούσες να καλύψεις τα σημεία με ένα τεντομένο σεντόνι. Τα simplex είναι μια οικογένεια τέτοιων σχημάτων, απο την οποία οι κανονικοί της εκπρόσωποι είναι το μοναδιαίο ευθυγραμμό τμήμα [0,1] στην 1 διάσταση,

και το τρίγωνο [(0,0), (0,1), (1,0)] στις 2 διαστάσεις

Κανονικοί εκπρόσωποι σε μεγαλύτερες διαστάσεις υπάρχουν, και συμβολίζονται με [0, e_1, e_2, ... e_n], όπου το 0 είναι n-διάστατο διάνυσμα, και τα e_i είναι παντού μηδέν εκτός από την i-οστή θέση που ειναι 1.

 

Ο λόγος που λέγονται κανονικοί εκπρόσωποι είναι επειδή οι κορυφές έχουν πολύ συγκεκριμένη θέση. Γενικά, τα 1-simplex είναι όλα τα ευθύγραμμα τμήματα, 2-simplex όλα τα τρίγωνα κλπ. Ο συμβολισμός παραμένει ο ίδιος: [v_1, v_2, ..., v_(n+1)] είναι το γενικό n-simplex.

 

Όμως οι κορυφές πρέπει να ικανοποιούν κάποια ιδιότητα ειδάλλως υπάρχουν περιπτώσεις εκφυλισμού, π.χ αν όλες οι κορυφές είναι σε μια γραμμή, τότε δεν παίρνεις τρίγωνο αλλά ευθύγραμμο τμήμα. Η ιδιότητα που εξασφαλίζει την δημιουργία του n-simplex ονομάζεται affine independence και μοιάζει με την linear independence (γραμμική ανεξαρτησία) (για αυτό και υλοποίησα την lindep() και όχι την affdep).

 

Η τελευταία παρατήρηση για τον υπολογισμό του όγκου είναι ότι 2 (πανομοιότυπα) τρίγωνα κάνουν 1 παραλληλόγραμμο. Η γενίκευση στις n διαστάσεις είναι ότι n! (παραγωντικό) simplices κάνουν ένα υπερπαραλληλεπίπεδο. Η διακρίνουσα υπολογίζει τον προσημασμένο όγκο του υπερπαραλληλεπιπέδου που σχηματίζει η συλλογή των διανυσμάτων. Επομένως, διαιρώντας με το παραγωντικό βρίσκουμε τον όγκο του simplex.

 

Π.χ. τα διανύσματα A,B,C φτιάχνουν το εξής παραλληλεπίπεδο

και, πίστεψε το ή δοκίμασε το, το παραπάνω παραλληλεπίπεδο χωράει ακριβώς 6 (= 3!) πανομοιότυπες τριγωνικές πυραμίδες (η κλειστή θήκη [0, A, B, C]... το σεντονάκι που λέγαμε)

 

V.I.Smirnov: Είσαι σίγουρος ότι αυτό είναι γρήγορο στην γενική περίπτωση (δηλ. d > 3); Αν και στο παραπάνω πρόγραμμα χρησιμοποιώ μια αργή διακρίνουσα, υπάρχουν γρήγοροι υπολογισμοί της διακρίνουσας και συνεπώς του όγκου του simplex. Ακούγεται ενδιαφέρον, έχεις κάνα λινκ (το graphic gems ίσως το έχει μόνο για d=3);

[V.I.Smirnov]

Τα simplex ("άπλοκα" στα ελληνικά) είναι η γενίκευση των τριγώνων σε περισσότερες από δυο διαστάσεις και

είναι μια ειδική κατηγορία κυρτών πολυγώνων.

Ειδικότερα, η convex hull ενός πεπερασμένου πλήθους σημείων στον R^d  λέγεται "πολύτοπο" (polytope).

Ένα πολύτοπο του οποίου η affine hull είναι διάστασης k λέγεται "k-πολύτοπο" (k-polytope).

( Affine hull είναι ένας υπόχωρος του R^d  όπου οι συντελεστές των γραμμικών συνδυασμών των

εν λόγω σημείων έχουν άθροισμα ακριβώς 1.)

H convex hull k+1 σημείων με k≤d και όχι σε affine space k-1, είναι ένα ειδικό k-πολύτοπο που

λέγεται simplex ή γενικότερα k-simplex.

Έτσι, σε δυο διαστάσεις, ένα 2-simplex είναι ένα τρίγωνο,

σε τρεις διαστάσεις, ένα 3-simplex είναι ένα τετράπλευρο.

Σε υψηλότερες διαστάσεις λέμε απλώς "simplex" ανεξάρτητα από τη διάσταση.

 

Τα παραπάνω τα θυμάμαι από ένα βιβλίο που αναφερόταν σε mesh generation.

Περισσότερες λεπτομέρειες υπάρχουν σε βιβλία τοπολογίας.

 

Αρκετές σχέσεις που αφορούν τρίγωνα και  τετράεδρα γενικεύονται εύκολα για simplex

(π.χ. μια ορίζουσα σε δυο διαστάσεις δίνει το εμβαδόν τριγώνου, σε 3 διαστάσεις τον όγκο του τετραέδρου κλπ),

γι αυτό έγραψα ότι η ιδέα δεν είναι εντυπωσιακή.

Σε 2 ή 3 ή διαστάσεις ο υπολογισμός των οριζουσών είναι εύκολος διότι υπάρχουν οι σχέσεις που τις δίνουν αναπτυγμένες.

Για μεγάλες ορίζουσες δυσκολεύουν τα πράγματα. Ξέρω ένα-δυο καλούς τρόπους αλλά δεν μπορώ να τους αναφέρω εδώ.

Στον ένα η ιδέα είναι να μετατραπεί ο αντίστοιχος πίνακας σε άνω τριγωνικό,

στον άλλο χρησιμοποιείται παρεμβολή.

 

Σε κοινές εφαρμογές γεωμετρίας χρειάζονται περισσότερο ο όγκος και η επιφάνεια τυχαίων πολυγώνων και πολυέδρων,

και αυτά υπολογίζονται σχετικά εύκολα με εφαρμογή των θ. Stokes και Gauss.

Στα graphics gems υπήρχαν τέτοιοι τύποι αλλά όπως είπες για 2 ή 3 διαστάσεις.

Τότε που τα διάβαζα (90'ς) ήμουν μαθητής και δεν ήξερα καλά τα θεωρήματα αυτά.

 

-

[filippos723]

υπαρχει κανενας που να εχει γνωσεις για κατασκευη hack σε online παιχνιδι ;;

καλη πληρωμη

[Kercyn]

Παιχνίδι με simplex; Ακούγεται πολύ βαρετό. Τι hack θέλεις, να προσθέτεις διαστάσεις;