Changeset 174 for cpp/frams/genetics


Ignore:
Timestamp:
03/14/14 00:45:03 (6 years ago)
Author:
oriona
Message:

Chambers in Foraminifera grow so that the length of the local communication path is minimal

Location:
cpp
Files:
2 added
4 edited

Legend:

Unmodified
Added
Removed
  • cpp

    • Property svn:ignore set to
      part_shapes
  • cpp/frams/genetics/fF/conv_fF.cpp

    r157 r174  
    88#include <common/nonstd_stl.h>
    99
    10 GenoConv_fF0::GenoConv_fF0()
    11 {
    12         name = "7-value Foraminifera encoding";
    13         in_format = 'F';
    14         out_format = '0';
    15         mapsupport = 0;
    16 }
    17 
    18 SString GenoConv_fF0::convert(SString &in, MultiMap *map)
    19 {
    20         fF_growth_params gp;
    21         if (!gp.load(in)) //invalid input genotype?
    22                 return ""; //so we return an invalid f0 genotype
    23         Model m;
    24         m.open();
    25         // subsequent parts (chambers) are placed relative to the previous part's orientation and location
    26         Part *p1, *p2;
    27         p1 = m.addNewPart(Part::SHAPE_ELLIPSOID); //initial part
    28         Pt3D scaling(gp.scalex, gp.scaley, gp.scalez);
    29         p1->scale = scaling;
    30 
    31         Orient rotation = Orient_1; //must be initialized explicitly because the default Orient constructor does not initialize anything
    32         rotation.rotate(Pt3D(0, gp.angle1, gp.angle2)); //assumed rotation around x is 0, which is limiting if we use assymetrical shapes (i.e., not spheres). But the original model had only two angles...
    33 
    34         for (int i = 1; i < gp.number_of_chambers; i++, p1 = p2)
    35         {
    36                 p2 = m.addNewPart(Part::SHAPE_ELLIPSOID);
    37                 p2->scale = p1->scale.entrywiseProduct(scaling); //each part's scale is its predecessor's scale * scaling
    38 
    39                 p2->setOrient(p1->o.transform(rotation)); //rotation transformed by p1 orientation
    40 
    41                 // Part's tip (offset from the center) is the point at (scale.x,0,0) in part's local coordinates.
    42                 // Get the offset in global coordinates by transforming it by the part's orientation.
    43                 p2->p = p1->p
    44                         + p1->o.transform(Pt3D(p1->scale.x*gp.translation,0,0))  // p1's tip transformed
    45                         + p2->o.transform(Pt3D(p2->scale.x*gp.translation,0,0)); // p2's tip transformed
    46 
    47                 m.addNewJoint(p1, p2, Joint::SHAPE_SOLID); //all parts must be connected
    48         }
    49         m.close();
    50         return m.getF0Geno().getGene();
    51 }
     10GenoConv_fF0::GenoConv_fF0() {
     11    name = "7-value Foraminifera encoding";
     12    in_format = 'F';
     13    out_format = '0';
     14    mapsupport = 0;
     15    cosines = new double[LATITUDE_NUM];
     16    sines = new double[LATITUDE_NUM]; 
     17    fill_cos_and_sin();
     18}
     19
     20GenoConv_fF0::~GenoConv_fF0() {
     21    delete []cosines;
     22    delete []sines;
     23}
     24
     25SString GenoConv_fF0::convert(SString &in, MultiMap *map) {
     26    fF_growth_params gp;
     27    if (!gp.load(in)) //invalid input genotype?
     28        return ""; //so we return an invalid f0 genotype
     29
     30    double div_radius_length = 1;//div_radius_lenght=1 or kx=ky=kz=1
     31    double radius = 1;
     32
     33    Model m;
     34    m.open();
     35    // subsequent parts (chambers) are placed relative to the previous part's orientation and location
     36    Part *p1, *p2;
     37
     38    fF_chamber3d **chambers = new fF_chamber3d*[gp.number_of_chambers];
     39
     40    for (int i = 0; i < gp.number_of_chambers; i++) {
     41        createSphere(i, chambers, radius, div_radius_length, gp.translation, gp.angle1, gp.angle2, gp.scalex, gp.scaley, gp.scalez);
     42    }
     43
     44    p1 = m.addNewPart(Part::SHAPE_ELLIPSOID);
     45    p1->p = Pt3D(chambers[0]->centerX, chambers[0]->centerY, chambers[0]->centerZ);
     46
     47
     48    for (int i = 1; i < gp.number_of_chambers; i++, p1 = p2) {
     49        p2 = m.addNewPart(Part::SHAPE_ELLIPSOID);
     50        p2->scale = p1->scale.entrywiseProduct(Pt3D(gp.scalex, gp.scaley, gp.scalez)); //each part's scale is its predecessor's scale * scaling
     51
     52        p2->p = Pt3D(chambers[i]->centerX, chambers[i]->centerY, chambers[i]->centerZ);
     53
     54        m.addNewJoint(p1, p2, Joint::SHAPE_SOLID); //all parts must be connected
     55    }
     56
     57    for (int i = 0; i < gp.number_of_chambers; i++)
     58        delete chambers[i];
     59    delete []chambers;
     60
     61    m.close();
     62    return m.getF0Geno().getGene();
     63}
     64
     65void GenoConv_fF0::createSphere(int which, fF_chamber3d **chambers, double radius_, double div_radius_length_, double div_vector_length_,
     66        double alpha_, double gamma_, double kx_, double ky_, double kz_)
     67{   
     68    chambers[which]=new fF_chamber3d(0.0f, 0.0f, 0.0f,
     69            (float) radius_, (float) radius_ * (float) kx_, 0.0f, 0.0f,
     70            (float) (radius_ * div_vector_length_), 0.0f, 0.0f, 0.0f, 0.0f);
     71    if (which == 0)
     72        chambers[which]->points = generate_points(chambers[which], which, kx_, ky_, kz_);
     73    if (which > 0) {
     74        /* old radius */
     75        double radiusOld, radius;
     76        radiusOld = chambers[which - 1]->radius;
     77        radius = div_radius_length_ * radiusOld;
     78        /* new growth vector length */
     79        double len = radius * div_vector_length_;
     80        if (radius < TOO_LITTLE) {
     81            radius = TOO_LITTLE;
     82            if (fabs(len) > (TOO_MUCH * radius)) {
     83                len = ((len < 0) ? (-1) : 1) * TOO_MUCH * radius;
     84            }
     85        }
     86        if (len == 0) {
     87            len = -0.0000001;
     88        }
     89
     90        /* aperture of the previous chamber */
     91        double pzx = chambers[which - 1]->holeX;
     92        double pzy = chambers[which - 1]->holeY;
     93        double pzz = chambers[which - 1]->holeZ;
     94
     95        //center of the previous chamber
     96        double pcx = chambers[which - 1]->centerX;
     97        double pcy = chambers[which - 1]->centerY;
     98        double pcz = chambers[which - 1]->centerZ;
     99
     100        /* aperture of the next to last chamber */
     101        double ppx;
     102        double ppy;
     103        double ppz;
     104
     105        if (which == 1) {
     106            ppx = pcx;
     107            ppy = pcy;
     108            ppz = pcz;
     109        } else {
     110            ppx = chambers[which - 2]->holeX;
     111            ppy = chambers[which - 2]->holeY;
     112            ppz = chambers[which - 2]->holeZ;
     113        }
     114
     115        double pzxprim = pzx - ppx;
     116        double pzyprim = pzy - ppy;
     117        double angle;
     118
     119        angle = atan2(pzyprim, pzxprim);
     120        double alpha = angle - alpha_;
     121
     122
     123        double gamma = chambers[which - 1]->phi + gamma_;
     124
     125        /* x */
     126        double wx = len * cos(alpha);
     127        /* y */
     128        double wy = len * sin(alpha);
     129        /* y */
     130        double wz = len * sin(alpha) * sin(gamma);
     131
     132        /*center of the new sphere*/
     133        double x = pzx + wx;
     134        double y = pzy + wy;
     135        double z = pzz + wz;
     136
     137        chambers[which]->centerX = (float) x;
     138        chambers[which]->centerY = (float) y;
     139        chambers[which]->centerZ = (float) z;
     140        chambers[which]->radius= (float) radius;
     141        chambers[which]->vectorTfX = (float) wx;
     142        chambers[which]->vectorTfY = (float) wy;
     143        chambers[which]->vectorTfZ = (float) wz;
     144        chambers[which]->beta = (float) alpha;
     145        chambers[which]->phi = (float) gamma;
     146
     147        chambers[which]->points = generate_points(chambers[which], which, kx_,  ky_,  kz_);
     148        search_hid(which, chambers, kx_,  ky_,  kz_);
     149        int pun;
     150        pun = find_hole(which, pzx, pzy, pzz, chambers, kx_,  ky_,  kz_);
     151
     152        chambers[which]->holeX = (float) chambers[which]->points[pun][0];
     153        chambers[which]->holeY = (float) chambers[which]->points[pun][1];
     154        chambers[which]->holeZ = (float) chambers[which]->points[pun][2];
     155    }
     156}
     157
     158void GenoConv_fF0::fill_cos_and_sin() {
     159    int i;
     160    double pi = acos(-1.0);
     161    double angle = pi / (((double) LATITUDE_NUM)*0.5);
     162    for (i = 0; i < LATITUDE_NUM; i++) {
     163        cosines[i] = cos((double) i * angle);
     164        sines[i] = sin((double) i * angle);
     165    }
     166}
     167
     168double** GenoConv_fF0::generate_points(fF_chamber3d *chamber, int which, double kx_, double ky_, double kz_) {
     169    float radius = chamber->radius;
     170    float cenx = chamber->centerX;
     171    float ceny = chamber->centerY;
     172    float cenz = chamber->centerZ;
     173
     174    double maxX = 0;
     175    double maxY = 0;
     176    double minX = 0;
     177    double minY = 0;
     178    double minZ = 0;
     179
     180    double kx = 1;
     181    double ky = 1;
     182    double kz = 1;
     183
     184    if (which > 0) {
     185        for (int kt = 1; kt < (which + 1); kt++) {
     186            kx = kx * kx_;
     187            ky = ky * ky_;
     188            kz = kz * kz_;
     189        }
     190    }
     191
     192    int i, j;
     193    double x, y, z;
     194   
     195    double **points = new double*[SIZE];
     196    for (int i = 0; i < SIZE; i++) {
     197        points[i] = new double[4];
     198    }
     199   
     200    for (i = 0; i < LONGITUDE_NUM; i++) {
     201        if (kx_ == 1 && ky_ == 1 && kz_ == 1) {
     202            y = ceny + radius * cosines[i];
     203        }
     204        else {
     205            y = ceny + ky * cosines[i];
     206        }
     207        for (j = 0; j < LATITUDE_NUM; j++) {
     208            if (kx_ == 1 && ky_ == 1 && kz_ == 1) {
     209                points[(i * LATITUDE_NUM) + j][0] = x = cenx + radius * cosines[j] * sines[i];
     210                points[(i * LATITUDE_NUM) + j][1] = y;
     211                points[(i * LATITUDE_NUM) + j][2] = z = cenz + radius * sines[j] * sines[i];
     212            } else {
     213                points[(i * LATITUDE_NUM) + j][0] = x = cenx + kx * cosines[j] * sines[i];
     214                points[(i * LATITUDE_NUM) + j][1] = y;
     215                points[(i * LATITUDE_NUM) + j][2] = z = cenz + kz * sines[j] * sines[i];
     216            }
     217
     218            points[(i * LATITUDE_NUM) + j][3] = 1.0;
     219            if (x < minX) minX = x;
     220            if (x > maxX) maxX = x;
     221            if (y < minY) minY = y;
     222            if (y > maxY) maxY = y;
     223
     224            if (z < minZ) minZ = z;
     225        };
     226    };
     227    return points;
     228
     229}
     230
     231double GenoConv_fF0::dist(double x1, double y1, double z1, double x2, double y2, double z2) {
     232    return sqrt((x2 - x1)*(x2 - x1) + (y2 - y1)*(y2 - y1) + (z2 - z1)*(z2 - z1));
     233}
     234
     235void GenoConv_fF0::search_hid(int nr, fF_chamber3d **spheres, double kx_, double ky_, double kz_) {
     236
     237    int i, j;
     238    if (nr != 0) {
     239        for (i = 0; i < nr; i++) {
     240            for (j = 0; j < AMOUNT; j++) {
     241                double X = spheres[nr]->points[j][0];
     242                double Y = spheres[nr]->points[j][1];
     243                double Z = spheres[nr]->points[j][2];
     244
     245                double srX0 = spheres[i]->centerX;
     246                double srY0 = spheres[i]->centerY;
     247                double srZ0 = spheres[i]->centerZ;
     248
     249                double a2;
     250                double b2;
     251                double c2;
     252
     253                if (kx_ != 1) {
     254                    a2 = (kx_ * kx_);
     255                } else {
     256                    a2 = (spheres[i]->radius * spheres[i]->radius);
     257                }
     258
     259                if (ky_ != 1) {
     260                    b2 = (ky_ * ky_);
     261
     262                } else {
     263                    b2 = (spheres[i]->radius * spheres[i]->radius);
     264                }
     265
     266                c2 = (kz_ * spheres[i]->radius) * (kz_ * spheres[i]->radius);
     267
     268                double up1 = (X - srX0) * (X - srX0);
     269                double up2 = (Y - srY0) * (Y - srY0);
     270                double up3 = (Z - srZ0) * (Z - srZ0);
     271
     272                double exp = up1 / a2;
     273                double exp2 = up2 / b2;
     274                double exp3 = up3 / c2;
     275
     276                double result = exp + exp2 + exp3;
     277
     278                if (result < (THICK_RATIO)
     279                        ) {
     280                    spheres[nr]->points[j][3] = 0;
     281                }
     282            }
     283        }
     284    }
     285}
     286
     287int GenoConv_fF0::find_hole(int which, double x, double y, double z, fF_chamber3d **chambers, double kx_, double ky_, double kz_) {
     288    int i;
     289    double distance;
     290    int found = 0;
     291    double dist_found = 0;
     292    int first = 1;
     293
     294    for (i = 0; i < AMOUNT; i++) {
     295        if (chambers[which]->points[i][3] != 0) //nie jest wewnatrz inne komory
     296        {
     297            distance = sqrt((chambers[which]->points[i][0] - x)*(chambers[which]->points[i][0] - x) +
     298                    (chambers[which]->points[i][1] - y)*(chambers[which]->points[i][1] - y)+
     299                    (chambers[which]->points[i][2] - z)*(chambers[which]->points[i][2] - z));
     300            if (first != 0) {
     301                found = i;
     302                dist_found = distance;
     303                first = 0;
     304            };
     305            if (distance < dist_found) {
     306                if (which != 0) {
     307                    bool good = true;
     308                    for (int j = 0; j < which; j++) {
     309                        {
     310                            double X = chambers[which]->points[i][0];
     311                            double Y = chambers[which]->points[i][1];
     312                            double Z = chambers[which]->points[i][2];
     313
     314                            double srX0 = chambers[j]->centerX;
     315                            double srY0 = chambers[j]->centerY;
     316                            double srZ0 = chambers[j]->centerZ;
     317
     318                            double a2 = (kx_ * chambers[j]->radius) * (kx_ * chambers[j]->radius);
     319                            double b2 = (ky_ * chambers[j]->radius) * (ky_ * chambers[j]->radius);
     320                            double c2 = (kz_ * chambers[j]->radius) * (kz_ * chambers[j]->radius);
     321
     322                            double up1 = (X - srX0) * (X - srX0);
     323                            double up2 = (Y - srY0) * (Y - srY0);
     324                            double up3 = (Z - srZ0) * (Z - srZ0);
     325
     326                            double exp1 = up1 / a2;
     327                            double exp2 = up2 / b2;
     328                            double exp3 = up3 / c2;
     329
     330                            double result = exp1 + exp2 + exp3;
     331                            if (result < 1.0)
     332                            {
     333                                good = false;
     334                            }
     335                        }
     336                    }
     337                    if (good) {
     338                        found = i;
     339                        dist_found = distance;
     340                    }
     341                }
     342            };
     343        };
     344    };
     345
     346    return (found);
     347}
  • cpp/frams/genetics/fF/conv_fF.h

    r140 r174  
    66#define _CONV_FF_H_
    77
     8#define TOO_MUCH 0.75
     9#define TOO_LITTLE 0.10
     10
     11#define HOLE_RADIUS 0.05f
     12#define LONGITUDE_NUM 69
     13
     14#define LATITUDE_NUM ((LONGITUDE_NUM - 1)*2)
     15#define AMOUNT ((LATITUDE_NUM)*(LONGITUDE_NUM))
     16
     17#define THICK_RATIO 0.95
     18
     19#define SIZE LONGITUDE_NUM * LATITUDE_NUM + LATITUDE_NUM
     20
    821#include <frams/util/multimap.h>
    922#include <frams/util/sstring.h>
    1023#include <frams/genetics/genoconv.h>
    11 
    12 
     24#include "fF_chamber3d.h"
    1325
    1426// The f9->f0 converter
    15 class GenoConv_fF0: public GenoConverter
    16 {
     27
     28class GenoConv_fF0 : public GenoConverter {
    1729public:
    18         GenoConv_fF0();
    19 
    20         //implementation of the GenoConverter method
    21         SString convert(SString &in, MultiMap *map);
     30    GenoConv_fF0();
     31    ~GenoConv_fF0();
     32    //implementation of the GenoConverter method
     33    SString convert(SString &in, MultiMap *map);
    2234
    2335protected:
    24         //auxiliary methods
    25         //...
     36    double* cosines;
     37    double* sines;
     38    void createSphere(int ktora, fF_chamber3d **chambers, double radius, double div_radius_length, double div_vector_length,
     39    double alpha, double gamma, double kx, double ky, double kz);
     40    double** generate_points(fF_chamber3d *chamber, int which, double kx, double ky, double kz);
     41    void fill_cos_and_sin();
     42    double dist(double x1, double y1, double z1, double x2, double y2, double z2);
     43    void search_hid(int nr, fF_chamber3d **spheres, double kx_, double ky_, double kz_);
     44    int find_hole(int which, double x, double y, double z, fF_chamber3d **chambers, double kx_, double ky_, double kz_);
    2645};
    2746
    2847#endif
     48
     49
  • cpp/frams/genetics/fF/fF_genotype.cpp

    r145 r174  
    1010        { "fF", 1, 7, "fF" },
    1111        { "n", 0, PARAM_CANOMITNAME, "number of chambers", "d 1 15 1", FIELD(number_of_chambers), },
    12         { "sx", 0, PARAM_CANOMITNAME, "scale x", "f 0.5 1.0 0.9", FIELD(scalex), },
    13         { "sy", 0, PARAM_CANOMITNAME, "scale y", "f 0.5 1.0 0.9", FIELD(scaley), },
    14         { "sz", 0, PARAM_CANOMITNAME, "scale z", "f 0.5 1.0 0.9", FIELD(scalez), },
    15         { "tr", 0, PARAM_CANOMITNAME, "translation factor", "f 0 1 0.5", FIELD(translation), }, //formally -1..1, but let's avoid redundancy and elliminate symmetrical constructs
     12        { "sx", 0, PARAM_CANOMITNAME, "scale x", "f 1.0 1.1 1.01", FIELD(scalex), },
     13        { "sy", 0, PARAM_CANOMITNAME, "scale y", "f 1.0 1.1 1.01", FIELD(scaley), },
     14        { "sz", 0, PARAM_CANOMITNAME, "scale z", "f 1.0 1.1 1.01", FIELD(scalez), },
     15        { "tr", 0, PARAM_CANOMITNAME, "translation factor", "f -1 1 0.5", FIELD(translation), },
    1616        { "a1", 0, PARAM_CANOMITNAME, "angle 1", "f -3.1415926 3.1415926 0", FIELD(angle1), },
    1717        { "a2", 0, PARAM_CANOMITNAME, "angle 2", "f -3.1415926 3.1415926 0", FIELD(angle2), },
Note: See TracChangeset for help on using the changeset viewer.