source: cpp/frams/genetics/fF/conv_fF.cpp @ 174

Last change on this file since 174 was 174, checked in by oriona, 9 years ago

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

  • Property svn:eol-style set to native
File size: 11.4 KB
Line 
1// This file is a part of the Framsticks GDK.
2// Copyright (C) 2002-2014  Maciej Komosinski and Szymon Ulatowski.  See LICENSE.txt for details.
3// Refer to http://www.framsticks.com/ for further information.
4
5#include "conv_fF.h"
6#include "fF_genotype.h"
7#include <frams/model/model.h>
8#include <common/nonstd_stl.h>
9
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}
Note: See TracBrowser for help on using the repository browser.