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

Last change on this file since 176 was 176, checked in by Maciej Komosinski, 6 years ago

Renamed #defines that conflicted with system-wide defines for Windows (added fF_ prefix), updated ranges for scaling parameters, updated the simplest genotype, formatted sources

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