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

Last change on this file since 256 was 256, checked in by Maciej Komosinski, 9 years ago

Chambers in Foraminifera now have the proper orientation in 3D; simplified source

  • Property svn:eol-style set to native
File size: 8.6 KB
Line 
1// This file is a part of the Framsticks GDK.
2// Copyright (C) 1999-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 <common/nonstd_stl.h>
8#include <common/Convert.h>
9
10GenoConv_fF0::GenoConv_fF0()
11{
12        name = "7-value 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)
39{
40        fF_growth_params gp;
41        if (!gp.load(in)) //invalid input genotype?
42                return ""; //so we return an invalid f0 genotype
43
44        double div_radius_length = 1; //div_radius_length=1 or kx=ky=kz=1
45        double radius = 1;
46
47        Model m;
48        m.open();
49
50        m.vis_style = "foram"; //dedicated visual look for Foraminifera
51
52        fF_chamber3d **chambers = new fF_chamber3d*[gp.number_of_chambers];
53        for (int i = 0; i < gp.number_of_chambers; i++)
54                createSphere(i, chambers, radius, div_radius_length, gp.translation, gp.angle1, gp.angle2, gp.scalex, gp.scaley, gp.scalez);
55
56        Part *p1 = addNewPart(&m, chambers[0]);
57        for (int i = 1; i < gp.number_of_chambers; i++)
58        {
59                Part *p2 = addNewPart(&m, chambers[i]);
60                p2->scale = p1->scale.entrywiseProduct(Pt3D(gp.scalex, gp.scaley, gp.scalez)); //each part's scale is its predecessor's scale * scaling
61                m.addNewJoint(p1, p2, Joint::SHAPE_SOLID); //all parts must be connected
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().getGene();
71}
72
73void GenoConv_fF0::createSphere(int which, fF_chamber3d **chambers, double radius_, double div_radius_length_, double div_vector_length_,
74        double alpha_, double gamma_, double kx_, double ky_, double kz_)
75{
76        chambers[which] = new fF_chamber3d(0.0f, 0.0f, 0.0f,
77                (float)radius_, (float)radius_ * (float)kx_, 0.0f, 0.0f,
78                (float)(radius_ * div_vector_length_), 0.0f, 0.0f, 0.0f, 0.0f);
79        if (which == 0)
80                chambers[which]->points = generate_points(chambers[which], which, kx_, ky_, kz_);
81        if (which > 0) {
82                /* old radius */
83                double radiusOld, radius;
84                radiusOld = chambers[which - 1]->radius;
85                radius = div_radius_length_ * radiusOld;
86                /* new growth vector length */
87                double len = radius * div_vector_length_;
88                if (radius < fF_TOO_LITTLE) {
89                        radius = fF_TOO_LITTLE;
90                        if (fabs(len) > (fF_TOO_MUCH * radius)) {
91                                len = ((len < 0) ? (-1) : 1) * fF_TOO_MUCH * radius;
92                        }
93                }
94                if (len == 0) {
95                        len = -0.0000001;
96                }
97
98                /* aperture of the previous chamber */
99                double pzx = chambers[which - 1]->holeX;
100                double pzy = chambers[which - 1]->holeY;
101                double pzz = chambers[which - 1]->holeZ;
102
103                //center of the previous chamber
104                double pcx = chambers[which - 1]->centerX;
105                double pcy = chambers[which - 1]->centerY;
106                double pcz = chambers[which - 1]->centerZ;
107
108                /* aperture of the next to last chamber */
109                double ppx;
110                double ppy;
111                double ppz;
112
113                if (which == 1)
114                {
115                        ppx = pcx;
116                        ppy = pcy;
117                        ppz = pcz;
118                }
119                else
120                {
121                        ppx = chambers[which - 2]->holeX;
122                        ppy = chambers[which - 2]->holeY;
123                        ppz = chambers[which - 2]->holeZ;
124                }
125
126                double pzxprim = pzx - ppx;
127                double pzyprim = pzy - ppy;
128                double angle;
129
130                angle = Convert::atan_2(pzyprim, pzxprim);
131                double alpha = angle - alpha_;
132
133
134                double gamma = chambers[which - 1]->phi + gamma_;
135
136                /* x */
137                double wx = len * cos(alpha);
138                /* y */
139                double wy = len * sin(alpha);
140                /* y */
141                double wz = len * sin(alpha) * sin(gamma);
142
143                /*center of the new sphere*/
144                double x = pzx + wx;
145                double y = pzy + wy;
146                double z = pzz + wz;
147
148                chambers[which]->centerX = (float)x;
149                chambers[which]->centerY = (float)y;
150                chambers[which]->centerZ = (float)z;
151                chambers[which]->radius = (float)radius;
152                chambers[which]->vectorTfX = (float)wx;
153                chambers[which]->vectorTfY = (float)wy;
154                chambers[which]->vectorTfZ = (float)wz;
155                chambers[which]->beta = (float)alpha;
156                chambers[which]->phi = (float)gamma;
157
158                chambers[which]->points = generate_points(chambers[which], which, kx_, ky_, kz_);
159                search_hid(which, chambers, kx_, ky_, kz_);
160                int pun;
161                pun = find_hole(which, pzx, pzy, pzz, chambers, kx_, ky_, kz_);
162
163                chambers[which]->holeX = (float)chambers[which]->points[pun].x;
164                chambers[which]->holeY = (float)chambers[which]->points[pun].y;
165                chambers[which]->holeZ = (float)chambers[which]->points[pun].z;
166        }
167}
168
169void GenoConv_fF0::precompute_cos_and_sin()
170{
171        int i;
172        double pi = acos(-1.0);
173        double angle = pi / (((double)fF_LATITUDE_NUM)*0.5);
174        for (i = 0; i < fF_LATITUDE_NUM; i++)
175        {
176                cosines[i] = cos((double)i * angle);
177                sines[i] = sin((double)i * angle);
178        }
179}
180
181fF_point* GenoConv_fF0::generate_points(fF_chamber3d *chamber, int which, double kx_, double ky_, double kz_)
182{
183        float radius = chamber->radius;
184        float cenx = chamber->centerX;
185        float ceny = chamber->centerY;
186        float cenz = chamber->centerZ;
187
188        double maxX = 0;
189        double maxY = 0;
190        double minX = 0;
191        double minY = 0;
192        double minZ = 0;
193
194        double kx = 1;
195        double ky = 1;
196        double kz = 1;
197
198        if (which > 0)
199        {
200                for (int kt = 1; kt < (which + 1); kt++)
201                {
202                        kx = kx * kx_;
203                        ky = ky * ky_;
204                        kz = kz * kz_;
205                }
206        }
207
208        bool all_k_ones = kx_ == 1 && ky_ == 1 && kz_ == 1;
209
210        double rx = all_k_ones ? radius : kx;
211        double ry = all_k_ones ? radius : ky;
212        double rz = all_k_ones ? radius : kz;
213
214        fF_point *points = new fF_point[fF_SIZE];
215
216        for (int i = 0; i < fF_LONGITUDE_NUM; i++)
217        {
218                double y = ceny + ry * cosines[i];
219
220                for (int j = 0; j < fF_LATITUDE_NUM; j++)
221                {
222                        double x = cenx + rx * cosines[j] * sines[i];
223                        double z = cenz + rz * sines[j] * sines[i];
224                        fF_point &p = points[(i * fF_LATITUDE_NUM) + j];
225                        p.x = x;
226                        p.y = y;
227                        p.z = z;
228                        p.inside = false;
229
230                        if (x < minX) minX = x;
231                        if (x > maxX) maxX = x;
232                        if (y < minY) minY = y;
233                        if (y > maxY) maxY = y;
234
235                        if (z < minZ) minZ = z;
236                }
237        }
238        return points;
239
240}
241
242double GenoConv_fF0::dist(double x1, double y1, double z1, double x2, double y2, double z2)
243{
244        return sqrt((x2 - x1)*(x2 - x1) + (y2 - y1)*(y2 - y1) + (z2 - z1)*(z2 - z1));
245}
246
247void GenoConv_fF0::search_hid(int nr, fF_chamber3d **spheres, double kx_, double ky_, double kz_)
248{
249        double kxsq = kx_*kx_;
250        double kysq = ky_*ky_;
251        double kzsq = kz_*kz_;
252
253        for (int i = 0; i < nr; i++)
254        {
255                double srX0 = spheres[i]->centerX;
256                double srY0 = spheres[i]->centerY;
257                double srZ0 = spheres[i]->centerZ;
258
259                double radsq = spheres[i]->radius * spheres[i]->radius;
260
261                double a2 = kx_ != 1 ? kxsq : radsq;
262                double b2 = ky_ != 1 ? kysq : radsq;
263                double c2 = kzsq * radsq;
264
265                for (int j = 0; j < fF_AMOUNT; j++)
266                {
267                        fF_point &p = spheres[nr]->points[j];
268
269                        double up1 = (p.x - srX0) * (p.x - srX0);
270                        double up2 = (p.y - srY0) * (p.y - srY0);
271                        double up3 = (p.z - srZ0) * (p.z - srZ0);
272
273                        double exp1 = up1 / a2;
274                        double exp2 = up2 / b2;
275                        double exp3 = up3 / c2;
276
277                        double result = exp1 + exp2 + exp3;
278
279                        if (result < fF_THICK_RATIO)
280                        {
281                                p.inside = true;
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{
289        int found = -1;
290        double distsq_found;
291
292        double kxsq = kx_*kx_;
293        double kysq = ky_*ky_;
294        double kzsq = kz_*kz_;
295
296        for (int i = 0; i < fF_AMOUNT; i++)
297        {
298                fF_point &p = chambers[which]->points[i];
299                if (!p.inside) //it is not inside another chamber
300                {
301                        double distancesq = (p.x - x)*(p.x - x) + (p.y - y)*(p.y - y) + (p.z - z)*(p.z - z);
302                        if (found < 0)
303                        {
304                                found = i;
305                                distsq_found = distancesq;
306                        }
307                        if (distancesq < distsq_found)
308                        {
309                                if (which != 0)
310                                {
311                                        bool good = true;
312                                        for (int j = 0; j < which && good; j++)
313                                        {
314                                                double srX0 = chambers[j]->centerX;
315                                                double srY0 = chambers[j]->centerY;
316                                                double srZ0 = chambers[j]->centerZ;
317
318                                                double radsq = chambers[j]->radius * chambers[j]->radius;
319
320                                                double a2 = kxsq * radsq;
321                                                double b2 = kysq * radsq;
322                                                double c2 = kzsq * radsq;
323
324                                                double up1 = (p.x - srX0) * (p.x - srX0);
325                                                double up2 = (p.y - srY0) * (p.y - srY0);
326                                                double up3 = (p.z - srZ0) * (p.z - srZ0);
327
328                                                double exp1 = up1 / a2;
329                                                double exp2 = up2 / b2;
330                                                double exp3 = up3 / c2;
331
332                                                double result = exp1 + exp2 + exp3;
333                                                if (result < 1.0)
334                                                {
335                                                        good = false;
336                                                }
337                                        }
338                                        if (good)
339                                        {
340                                                found = i;
341                                                distsq_found = distancesq;
342                                        }
343                                }
344                        }
345                }
346        }
347
348        return found;
349}
Note: See TracBrowser for help on using the repository browser.