source: cpp/frams/genetics/fF/fF_conv.cpp @ 978

Last change on this file since 978 was 779, checked in by Maciej Komosinski, 7 years ago

Unified file names of all files involved in genetic conversions and operations so that they start with "f<format>_"

  • Property svn:eol-style set to native
File size: 8.0 KB
Line 
1// This file is a part of Framsticks SDK.  http://www.framsticks.com/
2// Copyright (C) 1999-2018  Maciej Komosinski and Szymon Ulatowski.
3// See LICENSE.txt for details.
4
5#include "fF_conv.h"
6#include "fF_genotype.h"
7#include <common/nonstd_stl.h>
8#include <common/Convert.h>
9
10GenoConv_fF0::GenoConv_fF0()
11{
12        name = "10-parameter Foraminifera encoding";
13        in_format = 'F';
14        out_format = '0';
15        mapsupport = 0;
16        cosines = new double[fF_LATITUDE_NUM];
17        sines = new double[fF_LATITUDE_NUM];
18        precompute_cos_and_sin();
19}
20
21GenoConv_fF0::~GenoConv_fF0()
22{
23        delete[] cosines;
24        delete[] sines;
25}
26
27Part *GenoConv_fF0::addNewPart(Model *m, const fF_chamber3d* c)
28{
29        Part *part = m->addNewPart(Part::SHAPE_ELLIPSOID);
30        part->p = Pt3D(c->centerX, c->centerY, c->centerZ);
31        Pt3D hole = Pt3D(c->holeX, c->holeY, c->holeZ);
32        Orient o;
33        o.lookAt(part->p - hole);
34        part->setOrient(o);
35        return part;
36}
37
38SString GenoConv_fF0::convert(SString &in, MultiMap *map, bool using_checkpoints)
39{
40        fF_growth_params gp;
41        if (!gp.load(in.c_str())) //invalid input genotype?
42                return ""; //so we return an invalid f0 genotype
43
44        Model m;
45        m.open(using_checkpoints);
46
47        m.vis_style = "foram"; //dedicated visual look for Foraminifera
48
49        fF_chamber3d **chambers = new fF_chamber3d*[gp.number_of_chambers];
50        for (int i = 0; i < gp.number_of_chambers; i++)
51                createSphere(i, chambers, gp.radius0x, gp.radius0y, gp.radius0z, gp.translation, gp.angle1, gp.angle2, gp.scalex, gp.scaley, gp.scalez);
52
53        Part *p1 = addNewPart(&m, chambers[0]);
54        p1->scale = Pt3D(gp.radius0x, gp.radius0y, gp.radius0z); //size of the initial chamber
55        m.checkpoint();
56        for (int i = 1; i < gp.number_of_chambers; i++)
57        {
58                Part *p2 = addNewPart(&m, chambers[i]);
59                p2->scale = p1->scale.entrywiseProduct(Pt3D(gp.scalex, gp.scaley, gp.scalez)); //each part's scale is its predecessor's scale * scaling
60                m.addNewJoint(p1, p2, Joint::SHAPE_FIXED); //all parts must be connected
61                m.checkpoint();
62                p1 = p2;
63        }
64
65        for (int i = 0; i < gp.number_of_chambers; i++)
66                delete chambers[i];
67        delete[] chambers;
68
69        m.close();
70        return m.getF0Geno().getGenes();
71}
72
73void GenoConv_fF0::createSphere(int which, fF_chamber3d **chambers, double radius0x, double radius0y, double radius0z, double translation, double alpha_, double gamma_, double kx_, double ky_, double kz_)
74{
75        chambers[which] = new fF_chamber3d(0.0, 0.0, 0.0,
76                radius0x, radius0y, radius0z, radius0x * kx_, 0.0, 0.0,
77                radius0x * translation, 0.0, 0.0, 0.0, 0.0);
78        if (which == 0)
79                chambers[which]->points = generate_points(chambers[which]);
80        if (which > 0)
81        {
82                chambers[which]->radius_x = get_radius(chambers[which - 1]->radius_x, kx_, chambers[0]->radius_x);
83                chambers[which]->radius_y = get_radius(chambers[which - 1]->radius_y, ky_, chambers[0]->radius_y);
84                chambers[which]->radius_z = get_radius(chambers[which - 1]->radius_z, kz_, chambers[0]->radius_z);
85
86                /* new growth vector length */
87                double len = chambers[which]->radius_y * translation;
88                double max_radius = fF_TOO_MUCH * chambers[which]->radius_y;
89                if (fabs(len) > (max_radius))
90                        len = ((len < 0) ? (-1) : 1) * max_radius;
91                if (len == 0)
92                        len = -0.0000001;
93
94                /* aperture of the previous chamber */
95                double pzx = chambers[which - 1]->holeX;
96                double pzy = chambers[which - 1]->holeY;
97                double pzz = chambers[which - 1]->holeZ;
98
99                //center of the previous chamber
100                double pcx = chambers[which - 1]->centerX;
101                double pcy = chambers[which - 1]->centerY;
102                //double pcz = chambers[which - 1]->centerZ; //not used
103
104                /* aperture of the next to last chamber */
105                double ppx;
106                double ppy;
107                //double ppz; //not used
108
109                if (which == 1)
110                {
111                        ppx = pcx;
112                        ppy = pcy;
113                        //ppz = pcz;
114                }
115                else
116                {
117                        ppx = chambers[which - 2]->holeX;
118                        ppy = chambers[which - 2]->holeY;
119                        //ppz = chambers[which - 2]->holeZ;
120                }
121
122                double pzxprim = pzx - ppx;
123                double pzyprim = pzy - ppy;
124                double angle;
125
126                angle = Convert::atan_2(pzyprim, pzxprim);
127                double alpha = angle - alpha_;
128
129                double gamma = chambers[which - 1]->phi + gamma_;
130
131                double wx = len * cos(alpha);
132                double wy = len * sin(alpha);
133                double wz = len * sin(alpha) * sin(gamma);
134
135                /*center of the new sphere*/
136                double x = pzx + wx;
137                double y = pzy + wy;
138                double z = pzz + wz;
139
140                chambers[which]->centerX = x;
141                chambers[which]->centerY = y;
142                chambers[which]->centerZ = z;
143                chambers[which]->vectorTfX = wx;
144                chambers[which]->vectorTfY = wy;
145                chambers[which]->vectorTfZ = wz;
146                chambers[which]->beta = alpha;
147                chambers[which]->phi = gamma;
148
149                chambers[which]->points = generate_points(chambers[which]);
150                search_hid(which, chambers);
151                int pun;
152                pun = find_hole(which, pzx, pzy, pzz, chambers);
153                if (pun < 0) //should never happen
154                {
155                        logPrintf("GenoConv_fF0", "createSphere", LOG_ERROR, "find_hole(%d) returned %d", which, pun);
156                        pun = 0;
157                }
158                chambers[which]->holeX = chambers[which]->points[pun].x;
159                chambers[which]->holeY = chambers[which]->points[pun].y;
160                chambers[which]->holeZ = chambers[which]->points[pun].z;
161        }
162}
163
164double GenoConv_fF0::get_radius(double prev_radius, double scale, double start_radius)
165{
166        double radius = prev_radius * scale;
167        double min_radius = fF_TOO_LITTLE*start_radius;
168        if (radius < min_radius) {
169                radius = min_radius;
170        }
171
172        return radius;
173}
174
175void GenoConv_fF0::precompute_cos_and_sin()
176{
177        double angle = M_PI * 2 / fF_LATITUDE_NUM;
178        for (int i = 0; i < fF_LATITUDE_NUM; i++)
179        {
180                cosines[i] = cos(i * angle);
181                sines[i] = sin(i * angle);
182        }
183}
184
185fF_point* GenoConv_fF0::generate_points(fF_chamber3d *chamber)
186{
187        double cenx = chamber->centerX;
188        double ceny = chamber->centerY;
189        double cenz = chamber->centerZ;
190
191        double rx = chamber->radius_x;
192        double ry = chamber->radius_y;
193        double rz = chamber->radius_z;
194
195        fF_point *points = new fF_point[fF_SIZE];
196
197        for (int i = 0; i < fF_LONGITUDE_NUM; i++)
198        {
199                double y = ceny + ry * cosines[i];
200
201                for (int j = 0; j < fF_LATITUDE_NUM; j++)
202                {
203                        double x = cenx + rx * cosines[j] * sines[i];
204                        double z = cenz + rz * sines[j] * sines[i];
205                        fF_point &p = points[(i * fF_LATITUDE_NUM) + j];
206                        p.x = x;
207                        p.y = y;
208                        p.z = z;
209                        p.inside = false;
210                }
211        }
212        return points;
213
214}
215
216template<typename T> T Square(T x) { return x * x; }
217
218double GenoConv_fF0::dist(double x1, double y1, double z1, double x2, double y2, double z2)
219{
220        return sqrt(Square(x2 - x1) + Square(y2 - y1) + Square(z2 - z1));
221}
222
223void GenoConv_fF0::search_hid(int nr, fF_chamber3d **chambers)
224{
225        for (int i = 0; i < nr; i++)
226        {
227                fF_chamber3d *chamber = chambers[i];
228
229                double rx_sq = Square(chamber->radius_x);
230                double ry_sq = Square(chamber->radius_y);
231                double rz_sq = Square(chamber->radius_z);
232
233                for (int j = 0; j < fF_AMOUNT; j++)
234                {
235                        fF_point &p = chambers[nr]->points[j];
236
237                        double upx = Square(p.x - chamber->centerX);
238                        double upy = Square(p.y - chamber->centerY);
239                        double upz = Square(p.z - chamber->centerZ);
240
241                        double expx = upx / rx_sq;
242                        double expy = upy / ry_sq;
243                        double expz = upz / rz_sq;
244
245                        double result = expx + expy + expz;
246
247                        if (result < fF_THICK_RATIO)
248                        {
249                                p.inside = true;
250                        }
251                }
252        }
253}
254
255int GenoConv_fF0::find_hole(int which, double x, double y, double z, fF_chamber3d **chambers)
256{
257        int found = -1;
258        double distsq_found;
259
260        for (int i = 0; i < fF_AMOUNT; i++)
261        {
262                fF_point &p = chambers[which]->points[i];
263                if (!p.inside) //it is not inside another chamber
264                {
265                        double distancesq = Square(p.x - x) + Square(p.y - y) + Square(p.z - z);
266                        if (found < 0)
267                        {
268                                found = i;
269                                distsq_found = distancesq;
270                        }
271                        if (distancesq < distsq_found)
272                        {
273                                if (which != 0)
274                                {
275                                        bool good = true;
276                                        for (int j = 0; j < which && good; j++)
277                                        {
278                                                fF_chamber3d *chamber = chambers[j];
279
280                                                double rx_sq = Square(chamber->radius_x);
281                                                double ry_sq = Square(chamber->radius_y);
282                                                double rz_sq = Square(chamber->radius_z);
283
284                                                double upx = Square(p.x - chamber->centerX);
285                                                double upy = Square(p.y - chamber->centerY);
286                                                double upz = Square(p.z - chamber->centerZ);
287
288                                                double expx = upx / rx_sq;
289                                                double expy = upy / ry_sq;
290                                                double expz = upz / rz_sq;
291
292                                                double result = expx + expy + expz;
293                                                if (result < 1.0)
294                                                {
295                                                        good = false;
296                                                }
297                                        }
298                                        if (good)
299                                        {
300                                                found = i;
301                                                distsq_found = distancesq;
302                                        }
303                                }
304                        }
305                }
306        }
307
308        return found;
309}
Note: See TracBrowser for help on using the repository browser.