Tableau trop grand pour passer dans la RAM ? Traitement de voxels

Bonjour,

j'ai une petite question.

j'ai une erreur lorsque je fais tourner mon programme en c++. Il s'agit de segmentatiion fault (core dumped).
sur internet je vois que c' est quand on essaye d'acceder a un espace memoire deja alloue pour autre chose.

mais j'ai l'impression que dans mon cas c'est plus "il n'y a plus assez de memoire"
c'est possible ca ? parce que je vois rien a ce sujet sur internet.

parce que quand je fais un test avec une matrice test [200][200][200] ca passe.
mais quand je fais des test avec la matrice sur laquelle je veux travailler [220][1750][1300] ca ne passe plus.

est-ce que ce type d'erreur vient toujours d'une mauvaise allocation memoire ou ca peut venir aussi du fait que le pc n'a pas assez de ram ?

Oui, là c'est clairement ton OS qui dit à ton programme : holà doucement, j'ai pas tant de place disponible pour toi !

En effet, ta matrice contient 500.500.000 cases ! supposons que tu mettes un int (soit 4 octets) dans chacune des cases, tu viens de demander 2Go de ram !

merci pour la reponse.

ba oui elle fait 4 Go la matrice il me semble (je suis en float).
mais j'ai 16 Go de ram.

apres j'ai cette matrice la. et la matrice calculee apres reconstruction fait 1700*1700*1300.
j'ai les deux matrices de definies.

mais bon je vois pas trop comment je pourrais faire autrement ... j'ai besoin de toutes les donnees de la matrice 1 pour realiser chaque voxel de la matrice 2.

a part avoir plus de ram y a pas d'autre solution ?

à moins de compiler en 64bits y'a pas moyen d'avoir plus de 2Go de ram pour un programme (3 en bidouillant bien windows).

Une solution peut être d'utiliser des algorithmes qui peuvent travailler sur des parties de matrice plutôt que sur la matrice entière.

je suis sous linux (soit ubuntu soit suze entreprise) depend du pc sur lequel je travaille.
il me semblait que sous linux il n'y avait pas cette limite de 2Go.

(je suis sous code::blocks sous linux 64 bits).

par contre il me semble que c'est pas parce que j'ai linux en 64 bits que je compile en 64 bits c'est ca ?

j'avais lu un truc la dessus mais je sais plus ou.

comment je pourrais utiliser toute ma ram ?

le probleme c'est que pour calculer le pixel [i][j][k] de la matrice 2.
j'utilise des formules qui utilisent toutes les valeurs des pixels de la matrice 1. (a peu de chose pret). je vois pas trop comment redecouper mes matrices.
a la limite la matrice 2 je pourrais essayer de la decouper en deux. ca serait toujours ca.
mais dans ce cas comment je fais pour la reformer a la fin ?
faut que je fasse une sorte de free memory de tout le reste. plus je concatene les deux et apres seulement je les sauvegardes ?

chaque voxel nécessite une formule à 500 millions de paramètres? ça paraît étonnant! C'est quel type de calcul? Tu ne peux pas calculer par partie?

Tu as vraiment besoin d'une matrice/tenseur aussi gros (bien qu'elle rentre dans ta ram) ? Tu alloues comment ta structure : statique ou dynamique ?

Comment alloues-tu, pour commencer? Deux cas peuvent t'empêcher (ou te gêner fortement) d'utiliser toute ta RAM:

1. Allocation sur la pile (valable seulement pour les petites structures) => il faut allouer sur le heap pour pouvoir utiliser toute la RAM
2. Allocation contigüe: si tu demandes des blocs continus de 2Go ou plus, il y a de fortes chance pour que le système n'en trouve pas. Il vaut mieux utiliser des structures de données allouées par fragments dans ce cas

L'endroit où tu peux considérer disposer de toute la mémoire que tu veux c'est ton disque dur :).

La mémoire est construite sur ce modèle d'ailleurs :
registres
caches
ram
hdd

le cpu va travailler avec les registres, s'il a besoin il va chercher dans le cache qui s'il a besoin va chercher dans la ram, et si besoin encore plus dans le hdd. (on peut étendre ce modèle en ajoutant des serveurs et du réseau, etc.)

Citation:

---

Envoyé par svagrim

par contre il me semble que c'est pas parce que j'ai linux en 64 bits que je compile en 64 bits c'est ca ?

---

Il faut utiliser une toolchain 64bits et même sur un linux 64bits tu peux avoir une version 32 bits installée

Citation:

---

Envoyé par therwald

Comment alloues-tu, pour commencer? Deux cas peuvent t'empêcher (ou te gêner fortement) d'utiliser toute ta RAM:

1. Allocation sur la pile (valable seulement pour les petites structures) => il faut allouer sur le heap pour pouvoir utiliser toute la RAM
2. Allocation contigüe: si tu demandes des blocs continus de 2Go ou plus, il y a de fortes chance pour que le système n'en trouve pas. Il vaut mieux utiliser des structures de données allouées par fragments dans ce cas

---

Sur 16Go de ram, si il y a pas trop de truc qui tourne, j'espère que Linux arrive à faire son boulot correctement et avoir un bloc de 2Go en continue (en supposant qu'il gère de si gros bloc d'un coup, de mémoire je crois que Linux a pas de limite).

reconstruction d'imagerie haute resolution par rayon x.
on simule le trajet que prend le rayons x depuis une source cone-beam jusqu'au voxel. en prenant en compte les probabilites que le rayon x ait pu interagir facon scattered effect dans la matiere. c' est une methode iterative. donc on fait une premiere estimation de la matiere.
et on reevalue la quantite d'interaction qu'a pu avoir le photon en fonction de la matiere qu'il est cense avoir traversee etc ... on estime donc l'energie du spectre en fonction de la quantite de matiere traversee et du type de matiere etc pour faire une correction et passer a l'iteration suivante.

le probleme c'est que vu les angles des faisceaux qui sont calcules dans le programme, je vois vraiment pas comment decouper ma matrice.

pour l'initialisation de la matrice j'ai essaye un tableau en statique :

Code:

---

float matrice2 [dimensions];

---

et

Code:

---

vector<float> matrice2(dimensions,0);

---

qui avec cette methode ca marche a peut pret jusqu'au segmentation fault en cours de route.
la premiere methode j' ai le segmentation fault des le debut. (des l' initialisation)

je voulais voir si je pouvais pas utiliser une matrice sparse (parce que je reconstruit un cylindre donc j'ai quelques zeros. mais je comprend pas trop comment l'utiliser ...

Bonjour,

Ton rayon passe plan par plan et tu estimes les réfractions suivant l'angle de pénétration de ta surface. Un truc dans le genre si j'ai bien compris non?

Dans ce cas, pourquoi ne pas te contenter d'un simple tableau a 2 dimensions et prendre un petit tableau à coté dans lequel tu entres les valeurs de tes calculs. Ainsi tu gagnes en mémoire ;)
Enfin, je ne sais plus qui m'a dit ça un jour sur ce forum mais ça donnait:

"Quand tu n'as pas assez de place, revois ta façon de penser et donc de coder" :ccool:

Citation:

---

Envoyé par Flob90

Sur 16Go de ram, si il y a pas trop de truc qui tourne, j'espère que Linux arrive à faire son boulot correctement et avoir un bloc de 2Go en continue

---

Je vais modérer mon propos: le risque de ne pas trouver de bloc contigü est non négligeable: 2 fois 2Go sur 16 c'est quand même une fraction respectable.
Effectivement, s'il n'y a pas trop de trucs autour, ça doit passer.
@svagrim:
Oui, la notation float matrice[dimension] essaie d'allouer sur la pile, et là ça passe pas, il n'y a pas tant de place dans la pile. Le vector alloue sur le heap, donc là tu as une chance. Mais en sachant qu'un vector consomme plus que sa taille si on le fait grandir progressivement, vu qu'à chaque agrandissement il doit allouer le bloc de mémoire pour la nouvelle taille (plus grande) et recopier les données avant de libérer l'ancienne mémoire de stockage. Et là on se rapproche donc du risque de pénurie...
A priori la solution est de réserver la taille avant.

Salut,

Personnellement, je tenterais d'abord, et surtout, de voir s'il n'y a pas possibilité de se contenter d'un ou deux (allez, montons jusqu'à trois) plans successifs.

En effet, l'utilisation de tableaux à plusieurs dimensions fait que, plus tu augmentera le nombre de la dimension "la plus à gauche", plus l'espace mémoire dont tu auras besoin sera important, c'est mathématique (vu qu'il faut à chaque fois trouver un espace mémoire qui correspond à la multiplication des dimensions qui se trouvent à droite):D

si je propose cela, c'est simplement parce que j'imagine que si tu traverses 200 (1700 ou meme 3 000 000 de) plans, la valeur d'origine (comprends: avant de traverser le plan) n'est j'aimais que le résultat de la traversée du plan précédant.

Tu peux donc envisager (ou, plutôt, tu devrais donc pouvoir le faire :D) de traverser les différents plans exactement de la manière dont tu pelles un oignon: couche après couche, et tu pourrais mettre à profit le temps de la simulation d'un plan particulier pour charger (dans un autre thread, dédié à cela) le plan suivant à traiter.

Le fait d'avoir plusieurs plan n'étant là que pour "servir de tampon" afin de limiter au maximum les latences dues à la lecture des informations depuis le disque dur ;)

Ce que je te propose là n'est jamais qu'une adaptation fort simple du principe de "double (triple)" buffering si cher aux concepteur de jeux vidéo :D

Mais l'énorme avantage, c'est que dans l'hypothèse pessimiste (1700 plans de 1700*1300), tu n'as besoin que de 26 520 000 bytes de mémoire (en comptant 4 byte par points) qui peuvent même ne pas forcément être contigus en mémoire (en fait, j'aurais du parler de 3*8 840 000 bytes :D), ce qui reste parfaitement dans les limites de ce qu'un système, même 32 bits, est amplement en mesure de fournir à une application ;)

non les reconstructions tomographiques ne sont pas forcement 2D (plan par plan).
c'est meme rarement le cas de nos jours en dehors des acquisitions synchrotrons ou les faisceaux peuvent etre quasi-parallele.
la majorite des sources etant dites cone-beams.

on a donc une geometrie 3d qui est tres difficilement reduisible.

j'essaye de comprendre bien ce que tu dis koala et je te reponds :)
la c'est pas tres clair dans ma tete tes histoires de plans.
comment tu charges en meme temps dans un autre thread ? c'est de la programmation multi-processeurs dont tu me parles ou un truc comme ca ?

Citation:

---

Envoyé par svagrim

non les reconstructions tomographiques ne sont pas forcement 2D (plan par plan).
c'est meme rarement le cas de nos jours en dehors des acquisitions synchrotrons ou les faisceaux peuvent etre quasi-parallele.
la majorite des sources etant dites cone-beams.

on a donc une geometrie 3d qui est tres difficilement reduisible.

j'essaye de comprendre bien ce que tu dis koala et je te reponds :)
la c'est pas tres clair dans ma tete tes histoires de plans.

---

Bien, tu parlais d'un tableau à trois dimensions.

C'est donc que tu as trois axes sur lesquels te "balader"

1. l'axe des X (qui correspond à la dimension la plus à droite de ton tableau) pour l'horizontale
2. l'axe des Y (qui correspond à la dimension centrale de ton tableau) pour la "verticale" et
3. l'axe des Z (qui correspond à la dimension de gauche de ton tableau) pour la "profondeur"

En utilisant le nom de l'axe, ton tableau correspond donc à test[Z][Y][X](pour faire simple)

Chaque index donné pour Z correspond donc à un "plan" qui peut être envisagé de manière plus ou moins distincte.

Si tu as deux plans qui arrivent l'un derrière l'autre, tu peux parfaitement afficher / traiter / manipuler les données du premier pendant que tu charge celle du plan qui vient juste derrière.

Citation:

---

comment tu charges en meme temps dans un autre thread ? c'est de la programmation multi-processeurs dont tu me parles ou un truc comme ca ?

---

Je pensais au multi threading, à vrai dire.

L'idée est que ton rayon va, pour chaque plan qu'il traverse en un point donné, réagir d'une manière spécifique à la matière qu'il traverse et que cette réaction va influencer directement la manière dont il réagira lorsqu'il traversera le plan suivant, si j'ai bien compris :question:

L'idée est donc de travailler avec uniquement trois plans:

1. Un plan qui correspond à l'état du rayon avant de traverser le plan (mais peut etre après avoir traversé le plan précédant),
2. Un plan qui correspond à ce pour quoi tu dois évaluer la réaction de ton rayon et
3. Un plan dans lequel tu peux déjà commencer à placer les données qui correspondent au plan suivant, histoire de ne pas perdre trop de temps pour la suite.

Tu utilises l'état du rayon qui se trouve dans le plan (1) pour évaluer la réaction du rayon lorsqu'il traverse le plan (2), et, dans le même temps, tu charges un thread de placer dans le plan (3) les données qui correspondent au plan suivant.

Une fois que tu as traité toutes les données du plan (2) ET que tu as chargé toutes les données du plan (3), tu te retrouves donc dans le plan (1) avec l'état du rayon après avoir traversé le plan (2).

Tu peux donc commencer à évaluer l'état du rayon lorsqu'il traverse le plan (3) pendant que tu charges les données du plan rayon dans le plan (2). En fait, tu ne fais qu'intervertir les plans (2) et (3) pour savoir sur quel plan tu travailles et quel plan sert à prendre les données du plan suivant ;)

De cette manière, tu peux traiter énormément de plans différents (le rayon peut traverser 120 000 plans si tu le souhaites) en n'ayant jamais besoin que de deux plans pour travailler : un plan "en cours de traitement" et un autre "en cours de chargement".

Mais cette solution n'est applicable que si tu peux travailler "tranche par tranche" (ou plutôt plan par plan).

Une autre solution serait peut etre de mettre simplement chaque point en relation avec ceux qui l'entourent, sans doute sous une forme proche d'un octree.

La structure est, certes, plus imposante, mais chaque noeud peut être géré de manière séparée en mémoire, et cela facilitera sans doute le fait de trouver l'espace mémoire pour représenter tous les noeuds, vu que l'espace mémoire ne doit alors plus être contigu ;)

Citation:

---

Envoyé par svagrim

la majorite des sources etant dites cone-beams.

on a donc une geometrie 3d qui est tres difficilement reduisible.
?

---

Donc on est bien dans une problématique de géométrie? Pas de problèmes de diffration ou autres en jeu? Optique géométrique point barre? Dans ce cas ne pourrais-tu pas ramener tes couches de matières à des plans par transformation, même si dans l'espace réel ce sont plutôt des portions de sphères?
EDIT: je voulais dire pas de phénomènes ondulatoires pour compliquer le tableau?

Salut

J'ai pas tout lu, je commenterais pas les solutions techniques proposées. Tes remarques me font craindre un problème majeur de performance. Je me demande si ton programme sera fonctionnel, même après avoir résolu le problème de tableau de 4 Go.

1. tu dois charger la totalité des données en mémoire avant de commencer à travailler : perte de temps

2. si j'ai bien compris, tu dois parcourir toutes les correspondances entre 2 tableaux 3D. Donc un algo en O(N^6) ?
Si pour des tableaux de 400*400*400, ça te prend 1 heure. Alors pour des tableaux de 800*800*800, ça prendra 64 heures. Et pour un tableau de 1600*1600*1600, ça te prendra 4k heures (5 mois de calcul...)

3. les performances seront encore pire que ça, les données ne seront coalescentes, tu auras des problèmes de cache, plus de vectorisation

4. pas de parallélisation possible à priori avec ton algo et ta structurelle de données actuelle

Bref, à mon sens, il faut réussir dans tous les cas à découper le problème en plusieurs sous problème plus petit à gérer. Donc repartir du début : revenir au principe physique pour trouver des moyens de réduire une problématique globale à des problématique locale ? Changer la structure de données pour utiliser des algos de segmentation de l'espace ? Choisir un algorithme permettant de travailler sur des fragments de données ? (ce qui permettrait de paralléliser sur gpu ou cluster)

Je suppose que tu regardé la littérature ? Il y a des algos parallèles de reconstruction d'irm je pense

Bon courage

D'un point de vue algo, tu n'as pas besoin de définir tout ton espace de travail.

Si tu ne te focalises que sur la zone d'arrivée de ton rayon X et ensuite de sa propagation après réflexion tu peux te restreindre à un calcul d'angle et une direction (du ricochet), un calcul de probabilité pour le plan suivant et passer à l'étape d'après.

Je trouve que penser à représenter un rubiks cube entier alors que tu ne t'intéresses qu'à un seul petit carré n'est pas forcément très utile. Dans ton cas, c'est la même chose vu que ton rayon X est tellement localisé que tu peux pratiquement en faire une approche à quelques points (l'équivalent d'environ 10 fois la distance entre 2 atomes par exemple)

Sur ce, j'espère que tu pourras trouver une solution plus facile à mettre en oeuvre :ccool: