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

Last change on this file since 1181 was 1130, checked in by Maciej Komosinski, 3 years ago

Used std::min(), std::max() explicitly to avoid compiler confusion. Used std::size() explicitly instead of the equivalent macro

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