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

Last change on this file since 177 was 177, checked in by Maciej Komosinski, 8 years ago

Optimizations and simplifications

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