|
Blender
V2.59
|
00001 /* mball.c 00002 * 00003 * MetaBalls are created from a single Object (with a name without number in it), 00004 * here the DispList and BoundBox also is located. 00005 * All objects with the same name (but with a number in it) are added to this. 00006 * 00007 * texture coordinates are patched within the displist 00008 * 00009 * $Id: mball.c 36332 2011-04-26 07:17:21Z campbellbarton $ 00010 * 00011 * ***** BEGIN GPL LICENSE BLOCK ***** 00012 * 00013 * This program is free software; you can redistribute it and/or 00014 * modify it under the terms of the GNU General Public License 00015 * as published by the Free Software Foundation; either version 2 00016 * of the License, or (at your option) any later version. 00017 * 00018 * This program is distributed in the hope that it will be useful, 00019 * but WITHOUT ANY WARRANTY; without even the implied warranty of 00020 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00021 * GNU General Public License for more details. 00022 * 00023 * You should have received a copy of the GNU General Public License 00024 * along with this program; if not, write to the Free Software Foundation, 00025 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. 00026 * 00027 * The Original Code is Copyright (C) 2001-2002 by NaN Holding BV. 00028 * All rights reserved. 00029 * 00030 * Contributor(s): Jiri Hnidek <jiri.hnidek@vslib.cz>. 00031 * 00032 * ***** END GPL LICENSE BLOCK ***** 00033 */ 00034 00040 #include <stdio.h> 00041 #include <string.h> 00042 #include <math.h> 00043 #include <stdlib.h> 00044 #include <ctype.h> 00045 #include <float.h> 00046 #include "MEM_guardedalloc.h" 00047 00048 #include "DNA_material_types.h" 00049 #include "DNA_object_types.h" 00050 #include "DNA_meta_types.h" 00051 #include "DNA_scene_types.h" 00052 00053 00054 #include "BLI_blenlib.h" 00055 #include "BLI_math.h" 00056 #include "BLI_utildefines.h" 00057 00058 00059 00060 #include "BKE_global.h" 00061 #include "BKE_main.h" 00062 00063 /* #include "BKE_object.h" */ 00064 #include "BKE_animsys.h" 00065 #include "BKE_scene.h" 00066 #include "BKE_library.h" 00067 #include "BKE_displist.h" 00068 #include "BKE_mball.h" 00069 #include "BKE_object.h" 00070 #include "BKE_material.h" 00071 00072 /* Global variables */ 00073 00074 static float thresh= 0.6f; 00075 static int totelem=0; 00076 static MetaElem **mainb; 00077 static octal_tree *metaball_tree = NULL; 00078 /* Functions */ 00079 00080 void unlink_mball(MetaBall *mb) 00081 { 00082 int a; 00083 00084 for(a=0; a<mb->totcol; a++) { 00085 if(mb->mat[a]) mb->mat[a]->id.us--; 00086 mb->mat[a]= NULL; 00087 } 00088 } 00089 00090 00091 /* do not free mball itself */ 00092 void free_mball(MetaBall *mb) 00093 { 00094 unlink_mball(mb); 00095 00096 if(mb->adt) { 00097 BKE_free_animdata((ID *)mb); 00098 mb->adt = NULL; 00099 } 00100 if(mb->mat) MEM_freeN(mb->mat); 00101 if(mb->bb) MEM_freeN(mb->bb); 00102 BLI_freelistN(&mb->elems); 00103 if(mb->disp.first) freedisplist(&mb->disp); 00104 } 00105 00106 MetaBall *add_mball(const char *name) 00107 { 00108 MetaBall *mb; 00109 00110 mb= alloc_libblock(&G.main->mball, ID_MB, name); 00111 00112 mb->size[0]= mb->size[1]= mb->size[2]= 1.0; 00113 mb->texflag= MB_AUTOSPACE; 00114 00115 mb->wiresize= 0.4f; 00116 mb->rendersize= 0.2f; 00117 mb->thresh= 0.6f; 00118 00119 return mb; 00120 } 00121 00122 MetaBall *copy_mball(MetaBall *mb) 00123 { 00124 MetaBall *mbn; 00125 int a; 00126 00127 mbn= copy_libblock(mb); 00128 00129 BLI_duplicatelist(&mbn->elems, &mb->elems); 00130 00131 mbn->mat= MEM_dupallocN(mb->mat); 00132 for(a=0; a<mbn->totcol; a++) { 00133 id_us_plus((ID *)mbn->mat[a]); 00134 } 00135 mbn->bb= MEM_dupallocN(mb->bb); 00136 00137 mbn->editelems= NULL; 00138 mbn->lastelem= NULL; 00139 00140 return mbn; 00141 } 00142 00143 static void extern_local_mball(MetaBall *mb) 00144 { 00145 if(mb->mat) { 00146 extern_local_matarar(mb->mat, mb->totcol); 00147 } 00148 } 00149 00150 void make_local_mball(MetaBall *mb) 00151 { 00152 Main *bmain= G.main; 00153 Object *ob; 00154 int local=0, lib=0; 00155 00156 /* - only lib users: do nothing 00157 * - only local users: set flag 00158 * - mixed: make copy 00159 */ 00160 00161 if(mb->id.lib==NULL) return; 00162 if(mb->id.us==1) { 00163 mb->id.lib= NULL; 00164 mb->id.flag= LIB_LOCAL; 00165 new_id(&bmain->mball, (ID *)mb, NULL); 00166 extern_local_mball(mb); 00167 00168 return; 00169 } 00170 00171 for(ob= G.main->object.first; ob && ELEM(0, lib, local); ob= ob->id.next) { 00172 if(ob->data == mb) { 00173 if(ob->id.lib) lib= 1; 00174 else local= 1; 00175 } 00176 } 00177 00178 if(local && lib==0) { 00179 mb->id.lib= NULL; 00180 mb->id.flag= LIB_LOCAL; 00181 00182 new_id(&bmain->mball, (ID *)mb, NULL); 00183 extern_local_mball(mb); 00184 } 00185 else if(local && lib) { 00186 MetaBall *mbn= copy_mball(mb); 00187 mbn->id.us= 0; 00188 00189 for(ob= G.main->object.first; ob; ob= ob->id.next) { 00190 if(ob->data == mb) { 00191 if(ob->id.lib==NULL) { 00192 ob->data= mbn; 00193 mbn->id.us++; 00194 mb->id.us--; 00195 } 00196 } 00197 } 00198 } 00199 } 00200 00201 /* most simple meta-element adding function 00202 * don't do context manipulation here (rna uses) */ 00203 MetaElem *add_metaball_element(MetaBall *mb, const int type) 00204 { 00205 MetaElem *ml= MEM_callocN(sizeof(MetaElem), "metaelem"); 00206 00207 unit_qt(ml->quat); 00208 00209 ml->rad= 2.0; 00210 ml->s= 2.0; 00211 ml->flag= MB_SCALE_RAD; 00212 00213 switch(type) { 00214 case MB_BALL: 00215 ml->type = MB_BALL; 00216 ml->expx= ml->expy= ml->expz= 1.0; 00217 00218 break; 00219 case MB_TUBE: 00220 ml->type = MB_TUBE; 00221 ml->expx= ml->expy= ml->expz= 1.0; 00222 00223 break; 00224 case MB_PLANE: 00225 ml->type = MB_PLANE; 00226 ml->expx= ml->expy= ml->expz= 1.0; 00227 00228 break; 00229 case MB_ELIPSOID: 00230 ml->type = MB_ELIPSOID; 00231 ml->expx= 1.2f; 00232 ml->expy= 0.8f; 00233 ml->expz= 1.0; 00234 00235 break; 00236 case MB_CUBE: 00237 ml->type = MB_CUBE; 00238 ml->expx= ml->expy= ml->expz= 1.0; 00239 00240 break; 00241 default: 00242 break; 00243 } 00244 00245 BLI_addtail(&mb->elems, ml); 00246 00247 return ml; 00248 } 00255 void tex_space_mball(Object *ob) 00256 { 00257 DispList *dl; 00258 BoundBox *bb; 00259 float *data, min[3], max[3] /*, loc[3], size[3] */; 00260 int tot, doit=0; 00261 00262 if(ob->bb==NULL) ob->bb= MEM_callocN(sizeof(BoundBox), "mb boundbox"); 00263 bb= ob->bb; 00264 00265 /* Weird one, this. */ 00266 /* INIT_MINMAX(min, max); */ 00267 (min)[0]= (min)[1]= (min)[2]= 1.0e30f; 00268 (max)[0]= (max)[1]= (max)[2]= -1.0e30f; 00269 00270 dl= ob->disp.first; 00271 while(dl) { 00272 tot= dl->nr; 00273 if(tot) doit= 1; 00274 data= dl->verts; 00275 while(tot--) { 00276 /* Also weird... but longer. From utildefines. */ 00277 DO_MINMAX(data, min, max); 00278 data+= 3; 00279 } 00280 dl= dl->next; 00281 } 00282 00283 if(!doit) { 00284 min[0] = min[1] = min[2] = -1.0f; 00285 max[0] = max[1] = max[2] = 1.0f; 00286 } 00287 /* 00288 loc[0]= (min[0]+max[0])/2.0f; 00289 loc[1]= (min[1]+max[1])/2.0f; 00290 loc[2]= (min[2]+max[2])/2.0f; 00291 00292 size[0]= (max[0]-min[0])/2.0f; 00293 size[1]= (max[1]-min[1])/2.0f; 00294 size[2]= (max[2]-min[2])/2.0f; 00295 */ 00296 boundbox_set_from_min_max(bb, min, max); 00297 } 00298 00299 float *make_orco_mball(Object *ob, ListBase *dispbase) 00300 { 00301 BoundBox *bb; 00302 DispList *dl; 00303 float *data, *orco, *orcodata; 00304 float loc[3], size[3]; 00305 int a; 00306 00307 /* restore size and loc */ 00308 bb= ob->bb; 00309 loc[0]= (bb->vec[0][0]+bb->vec[4][0])/2.0f; 00310 size[0]= bb->vec[4][0]-loc[0]; 00311 loc[1]= (bb->vec[0][1]+bb->vec[2][1])/2.0f; 00312 size[1]= bb->vec[2][1]-loc[1]; 00313 loc[2]= (bb->vec[0][2]+bb->vec[1][2])/2.0f; 00314 size[2]= bb->vec[1][2]-loc[2]; 00315 00316 dl= dispbase->first; 00317 orcodata= MEM_mallocN(sizeof(float)*3*dl->nr, "MballOrco"); 00318 00319 data= dl->verts; 00320 orco= orcodata; 00321 a= dl->nr; 00322 while(a--) { 00323 orco[0]= (data[0]-loc[0])/size[0]; 00324 orco[1]= (data[1]-loc[1])/size[1]; 00325 orco[2]= (data[2]-loc[2])/size[2]; 00326 00327 data+= 3; 00328 orco+= 3; 00329 } 00330 00331 return orcodata; 00332 } 00333 00334 /* Note on mball basis stuff 2.5x (this is a can of worms) 00335 * This really needs a rewrite/refactor its totally broken in anything other then basic cases 00336 * Multiple Scenes + Set Scenes & mixing mball basis SHOULD work but fails to update the depsgraph on rename 00337 * and linking into scenes or removal of basis mball. so take care when changing this code. 00338 * 00339 * Main idiot thing here is that the system returns find_basis_mball() objects which fail a is_basis_mball() test. 00340 * 00341 * Not only that but the depsgraph and their areas depend on this behavior!, so making small fixes here isn't worth it. 00342 * - Campbell 00343 */ 00344 00345 00351 int is_basis_mball(Object *ob) 00352 { 00353 int len; 00354 00355 /* just a quick test */ 00356 len= strlen(ob->id.name); 00357 if( isdigit(ob->id.name[len-1]) ) return 0; 00358 return 1; 00359 } 00360 00361 /* return nonzero if ob1 is a basis mball for ob */ 00362 int is_mball_basis_for(Object *ob1, Object *ob2) 00363 { 00364 int basis1nr, basis2nr; 00365 char basis1name[32], basis2name[32]; 00366 00367 BLI_split_name_num(basis1name, &basis1nr, ob1->id.name+2, '.'); 00368 BLI_split_name_num(basis2name, &basis2nr, ob2->id.name+2, '.'); 00369 00370 if(!strcmp(basis1name, basis2name)) return is_basis_mball(ob1); 00371 else return 0; 00372 } 00373 00374 /* \brief copy some properties from object to other metaball object with same base name 00375 * 00376 * When some properties (wiresize, threshold, update flags) of metaball are changed, then this properties 00377 * are copied to all metaballs in same "group" (metaballs with same base name: MBall, 00378 * MBall.001, MBall.002, etc). The most important is to copy properties to the base metaball, 00379 * because this metaball influence polygonisation of metaballs. */ 00380 void copy_mball_properties(Scene *scene, Object *active_object) 00381 { 00382 Scene *sce_iter= scene; 00383 Base *base; 00384 Object *ob; 00385 MetaBall *active_mball = (MetaBall*)active_object->data; 00386 int basisnr, obnr; 00387 char basisname[32], obname[32]; 00388 00389 BLI_split_name_num(basisname, &basisnr, active_object->id.name+2, '.'); 00390 00391 /* XXX recursion check, see scene.c, just too simple code this next_object() */ 00392 if(F_ERROR==next_object(&sce_iter, 0, NULL, NULL)) 00393 return; 00394 00395 while(next_object(&sce_iter, 1, &base, &ob)) { 00396 if (ob->type==OB_MBALL) { 00397 if(ob!=active_object){ 00398 BLI_split_name_num(obname, &obnr, ob->id.name+2, '.'); 00399 00400 /* Object ob has to be in same "group" ... it means, that it has to have 00401 * same base of its name */ 00402 if(strcmp(obname, basisname)==0){ 00403 MetaBall *mb= ob->data; 00404 00405 /* Copy properties from selected/edited metaball */ 00406 mb->wiresize= active_mball->wiresize; 00407 mb->rendersize= active_mball->rendersize; 00408 mb->thresh= active_mball->thresh; 00409 mb->flag= active_mball->flag; 00410 } 00411 } 00412 } 00413 } 00414 } 00415 00425 Object *find_basis_mball(Scene *scene, Object *basis) 00426 { 00427 Scene *sce_iter= scene; 00428 Base *base; 00429 Object *ob,*bob= basis; 00430 MetaElem *ml=NULL; 00431 int basisnr, obnr; 00432 char basisname[32], obname[32]; 00433 00434 BLI_split_name_num(basisname, &basisnr, basis->id.name+2, '.'); 00435 totelem= 0; 00436 00437 /* XXX recursion check, see scene.c, just too simple code this next_object() */ 00438 if(F_ERROR==next_object(&sce_iter, 0, NULL, NULL)) 00439 return NULL; 00440 00441 while(next_object(&sce_iter, 1, &base, &ob)) { 00442 00443 if (ob->type==OB_MBALL) { 00444 if(ob==bob){ 00445 MetaBall *mb= ob->data; 00446 00447 /* if bob object is in edit mode, then dynamic list of all MetaElems 00448 * is stored in editelems */ 00449 if(mb->editelems) ml= mb->editelems->first; 00450 /* if bob object is in object mode */ 00451 else ml= mb->elems.first; 00452 } 00453 else{ 00454 BLI_split_name_num(obname, &obnr, ob->id.name+2, '.'); 00455 00456 /* object ob has to be in same "group" ... it means, that it has to have 00457 * same base of its name */ 00458 if(strcmp(obname, basisname)==0){ 00459 MetaBall *mb= ob->data; 00460 00461 /* if object is in edit mode, then dynamic list of all MetaElems 00462 * is stored in editelems */ 00463 if(mb->editelems) ml= mb->editelems->first; 00464 /* if bob object is in object mode */ 00465 else ml= mb->elems.first; 00466 00467 if(obnr<basisnr){ 00468 if(!(ob->flag & OB_FROMDUPLI)){ 00469 basis= ob; 00470 basisnr= obnr; 00471 } 00472 } 00473 } 00474 } 00475 00476 while(ml){ 00477 if(!(ml->flag & MB_HIDE)) totelem++; 00478 ml= ml->next; 00479 } 00480 } 00481 } 00482 00483 return basis; 00484 } 00485 00486 00487 /* ******************** ARITH ************************* */ 00488 00489 /* BASED AT CODE (but mostly rewritten) : 00490 * C code from the article 00491 * "An Implicit Surface Polygonizer" 00492 * by Jules Bloomenthal, jbloom@beauty.gmu.edu 00493 * in "Graphics Gems IV", Academic Press, 1994 00494 00495 * Authored by Jules Bloomenthal, Xerox PARC. 00496 * Copyright (c) Xerox Corporation, 1991. All rights reserved. 00497 * Permission is granted to reproduce, use and distribute this code for 00498 * any and all purposes, provided that this notice appears in all copies. */ 00499 00500 #define RES 12 /* # converge iterations */ 00501 00502 #define L 0 /* left direction: -x, -i */ 00503 #define R 1 /* right direction: +x, +i */ 00504 #define B 2 /* bottom direction: -y, -j */ 00505 #define T 3 /* top direction: +y, +j */ 00506 #define N 4 /* near direction: -z, -k */ 00507 #define F 5 /* far direction: +z, +k */ 00508 #define LBN 0 /* left bottom near corner */ 00509 #define LBF 1 /* left bottom far corner */ 00510 #define LTN 2 /* left top near corner */ 00511 #define LTF 3 /* left top far corner */ 00512 #define RBN 4 /* right bottom near corner */ 00513 #define RBF 5 /* right bottom far corner */ 00514 #define RTN 6 /* right top near corner */ 00515 #define RTF 7 /* right top far corner */ 00516 00517 /* the LBN corner of cube (i, j, k), corresponds with location 00518 * (i-0.5)*size, (j-0.5)*size, (k-0.5)*size) */ 00519 00520 #define HASHBIT (5) 00521 #define HASHSIZE (size_t)(1<<(3*HASHBIT)) 00523 #define HASH(i,j,k) ((((( (i) & 31)<<5) | ( (j) & 31))<<5 ) | ( (k) & 31) ) 00524 00525 #define MB_BIT(i, bit) (((i)>>(bit))&1) 00526 #define FLIP(i,bit) ((i)^1<<(bit)) /* flip the given bit of i */ 00527 00528 00529 /* **************** POLYGONIZATION ************************ */ 00530 00531 void calc_mballco(MetaElem *ml, float *vec) 00532 { 00533 if(ml->mat) { 00534 mul_m4_v3((float ( * )[4])ml->mat, vec); 00535 } 00536 } 00537 00538 float densfunc(MetaElem *ball, float x, float y, float z) 00539 { 00540 float dist2 = 0.0, dx, dy, dz; 00541 float vec[3]; 00542 00543 vec[0]= x; 00544 vec[1]= y; 00545 vec[2]= z; 00546 mul_m4_v3((float ( * )[4])ball->imat, vec); 00547 dx= vec[0]; 00548 dy= vec[1]; 00549 dz= vec[2]; 00550 00551 if(ball->type==MB_BALL) { 00552 } 00553 else if(ball->type==MB_TUBEX) { 00554 if( dx > ball->len) dx-= ball->len; 00555 else if(dx< -ball->len) dx+= ball->len; 00556 else dx= 0.0; 00557 } 00558 else if(ball->type==MB_TUBEY) { 00559 if( dy > ball->len) dy-= ball->len; 00560 else if(dy< -ball->len) dy+= ball->len; 00561 else dy= 0.0; 00562 } 00563 else if(ball->type==MB_TUBEZ) { 00564 if( dz > ball->len) dz-= ball->len; 00565 else if(dz< -ball->len) dz+= ball->len; 00566 else dz= 0.0; 00567 } 00568 else if(ball->type==MB_TUBE) { 00569 if( dx > ball->expx) dx-= ball->expx; 00570 else if(dx< -ball->expx) dx+= ball->expx; 00571 else dx= 0.0; 00572 } 00573 else if(ball->type==MB_PLANE) { 00574 if( dx > ball->expx) dx-= ball->expx; 00575 else if(dx< -ball->expx) dx+= ball->expx; 00576 else dx= 0.0; 00577 if( dy > ball->expy) dy-= ball->expy; 00578 else if(dy< -ball->expy) dy+= ball->expy; 00579 else dy= 0.0; 00580 } 00581 else if(ball->type==MB_ELIPSOID) { 00582 dx *= 1/ball->expx; 00583 dy *= 1/ball->expy; 00584 dz *= 1/ball->expz; 00585 } 00586 else if(ball->type==MB_CUBE) { 00587 if( dx > ball->expx) dx-= ball->expx; 00588 else if(dx< -ball->expx) dx+= ball->expx; 00589 else dx= 0.0; 00590 if( dy > ball->expy) dy-= ball->expy; 00591 else if(dy< -ball->expy) dy+= ball->expy; 00592 else dy= 0.0; 00593 if( dz > ball->expz) dz-= ball->expz; 00594 else if(dz< -ball->expz) dz+= ball->expz; 00595 else dz= 0.0; 00596 } 00597 00598 dist2= (dx*dx + dy*dy + dz*dz); 00599 00600 if(ball->flag & MB_NEGATIVE) { 00601 dist2= 1.0f-(dist2/ball->rad2); 00602 if(dist2 < 0.0f) return 0.5f; 00603 00604 return 0.5f-ball->s*dist2*dist2*dist2; 00605 } 00606 else { 00607 dist2= 1.0f-(dist2/ball->rad2); 00608 if(dist2 < 0.0f) return -0.5f; 00609 00610 return ball->s*dist2*dist2*dist2 -0.5f; 00611 } 00612 } 00613 00614 octal_node* find_metaball_octal_node(octal_node *node, float x, float y, float z, short depth) 00615 { 00616 if(!depth) return node; 00617 00618 if(z < node->z){ 00619 if(y < node->y){ 00620 if(x < node->x){ 00621 if(node->nodes[0]) 00622 return find_metaball_octal_node(node->nodes[0],x,y,z,depth--); 00623 else 00624 return node; 00625 } 00626 else{ 00627 if(node->nodes[1]) 00628 return find_metaball_octal_node(node->nodes[1],x,y,z,depth--); 00629 else 00630 return node; 00631 } 00632 } 00633 else{ 00634 if(x < node->x){ 00635 if(node->nodes[3]) 00636 return find_metaball_octal_node(node->nodes[3],x,y,z,depth--); 00637 else 00638 return node; 00639 } 00640 else{ 00641 if(node->nodes[2]) 00642 return find_metaball_octal_node(node->nodes[2],x,y,z,depth--); 00643 else 00644 return node; 00645 } 00646 } 00647 } 00648 else{ 00649 if(y < node->y){ 00650 if(x < node->x){ 00651 if(node->nodes[4]) 00652 return find_metaball_octal_node(node->nodes[4],x,y,z,depth--); 00653 else 00654 return node; 00655 } 00656 else{ 00657 if(node->nodes[5]) 00658 return find_metaball_octal_node(node->nodes[5],x,y,z,depth--); 00659 else 00660 return node; 00661 } 00662 } 00663 else{ 00664 if(x < node->x){ 00665 if(node->nodes[7]) 00666 return find_metaball_octal_node(node->nodes[7],x,y,z,depth--); 00667 else 00668 return node; 00669 } 00670 else{ 00671 if(node->nodes[6]) 00672 return find_metaball_octal_node(node->nodes[6],x,y,z,depth--); 00673 else 00674 return node; 00675 } 00676 } 00677 } 00678 00679 return node; 00680 } 00681 00682 float metaball(float x, float y, float z) 00683 /* float x, y, z; */ 00684 { 00685 struct octal_node *node; 00686 struct ml_pointer *ml_p; 00687 float dens=0; 00688 int a; 00689 00690 if(totelem > 1){ 00691 node= find_metaball_octal_node(metaball_tree->first, x, y, z, metaball_tree->depth); 00692 if(node){ 00693 ml_p= node->elems.first; 00694 00695 while(ml_p){ 00696 dens+=densfunc(ml_p->ml, x, y, z); 00697 ml_p= ml_p->next; 00698 } 00699 00700 dens+= -0.5f*(metaball_tree->pos - node->pos); 00701 dens+= 0.5f*(metaball_tree->neg - node->neg); 00702 } 00703 else{ 00704 for(a=0; a<totelem; a++) { 00705 dens+= densfunc( mainb[a], x, y, z); 00706 } 00707 } 00708 } 00709 else{ 00710 dens+= densfunc( mainb[0], x, y, z); 00711 } 00712 00713 return thresh - dens; 00714 } 00715 00716 /* ******************************************** */ 00717 00718 static int *indices=NULL; 00719 static int totindex, curindex; 00720 00721 00722 void accum_mballfaces(int i1, int i2, int i3, int i4) 00723 { 00724 int *newi, *cur; 00725 /* static int i=0; I would like to delete altogether, but I don't dare to, yet */ 00726 00727 if(totindex==curindex) { 00728 totindex+= 256; 00729 newi= MEM_mallocN(4*sizeof(int)*totindex, "vertindex"); 00730 00731 if(indices) { 00732 memcpy(newi, indices, 4*sizeof(int)*(totindex-256)); 00733 MEM_freeN(indices); 00734 } 00735 indices= newi; 00736 } 00737 00738 cur= indices+4*curindex; 00739 00740 /* displists now support array drawing, we treat tri's as fake quad */ 00741 00742 cur[0]= i1; 00743 cur[1]= i2; 00744 cur[2]= i3; 00745 if(i4==0) 00746 cur[3]= i3; 00747 else 00748 cur[3]= i4; 00749 00750 curindex++; 00751 00752 } 00753 00754 /* ******************* MEMORY MANAGEMENT *********************** */ 00755 void *new_pgn_element(int size) 00756 { 00757 /* during polygonize 1000s of elements are allocated 00758 * and never freed in between. Freeing only done at the end. 00759 */ 00760 int blocksize= 16384; 00761 static int offs= 0; /* the current free address */ 00762 static struct pgn_elements *cur= NULL; 00763 static ListBase lb= {NULL, NULL}; 00764 void *adr; 00765 00766 if(size>10000 || size==0) { 00767 printf("incorrect use of new_pgn_element\n"); 00768 } 00769 else if(size== -1) { 00770 cur= lb.first; 00771 while(cur) { 00772 MEM_freeN(cur->data); 00773 cur= cur->next; 00774 } 00775 BLI_freelistN(&lb); 00776 00777 return NULL; 00778 } 00779 00780 size= 4*( (size+3)/4 ); 00781 00782 if(cur) { 00783 if(size+offs < blocksize) { 00784 adr= (void *) (cur->data+offs); 00785 offs+= size; 00786 return adr; 00787 } 00788 } 00789 00790 cur= MEM_callocN( sizeof(struct pgn_elements), "newpgn"); 00791 cur->data= MEM_callocN(blocksize, "newpgn"); 00792 BLI_addtail(&lb, cur); 00793 00794 offs= size; 00795 return cur->data; 00796 } 00797 00798 void freepolygonize(PROCESS *p) 00799 { 00800 MEM_freeN(p->corners); 00801 MEM_freeN(p->edges); 00802 MEM_freeN(p->centers); 00803 00804 new_pgn_element(-1); 00805 00806 if(p->vertices.ptr) MEM_freeN(p->vertices.ptr); 00807 } 00808 00809 /**** Cubical Polygonization (optional) ****/ 00810 00811 #define LB 0 /* left bottom edge */ 00812 #define LT 1 /* left top edge */ 00813 #define LN 2 /* left near edge */ 00814 #define LF 3 /* left far edge */ 00815 #define RB 4 /* right bottom edge */ 00816 #define RT 5 /* right top edge */ 00817 #define RN 6 /* right near edge */ 00818 #define RF 7 /* right far edge */ 00819 #define BN 8 /* bottom near edge */ 00820 #define BF 9 /* bottom far edge */ 00821 #define TN 10 /* top near edge */ 00822 #define TF 11 /* top far edge */ 00823 00824 static INTLISTS *cubetable[256]; 00825 00826 /* edge: LB, LT, LN, LF, RB, RT, RN, RF, BN, BF, TN, TF */ 00827 static int corner1[12] = { 00828 LBN,LTN,LBN,LBF,RBN,RTN,RBN,RBF,LBN,LBF,LTN,LTF}; 00829 static int corner2[12] = { 00830 LBF,LTF,LTN,LTF,RBF,RTF,RTN,RTF,RBN,RBF,RTN,RTF}; 00831 static int leftface[12] = { 00832 B, L, L, F, R, T, N, R, N, B, T, F}; 00833 /* face on left when going corner1 to corner2 */ 00834 static int rightface[12] = { 00835 L, T, N, L, B, R, R, F, B, F, N, T}; 00836 /* face on right when going corner1 to corner2 */ 00837 00838 00839 /* docube: triangulate the cube directly, without decomposition */ 00840 00841 void docube(CUBE *cube, PROCESS *p, MetaBall *mb) 00842 { 00843 INTLISTS *polys; 00844 CORNER *c1, *c2; 00845 int i, index = 0, count, indexar[8]; 00846 00847 for (i = 0; i < 8; i++) if (cube->corners[i]->value > 0.0f) index += (1<<i); 00848 00849 for (polys = cubetable[index]; polys; polys = polys->next) { 00850 INTLIST *edges; 00851 00852 count = 0; 00853 00854 for (edges = polys->list; edges; edges = edges->next) { 00855 c1 = cube->corners[corner1[edges->i]]; 00856 c2 = cube->corners[corner2[edges->i]]; 00857 00858 indexar[count] = vertid(c1, c2, p, mb); 00859 count++; 00860 } 00861 if(count>2) { 00862 switch(count) { 00863 case 3: 00864 accum_mballfaces(indexar[2], indexar[1], indexar[0], 0); 00865 break; 00866 case 4: 00867 if(indexar[0]==0) accum_mballfaces(indexar[0], indexar[3], indexar[2], indexar[1]); 00868 else accum_mballfaces(indexar[3], indexar[2], indexar[1], indexar[0]); 00869 break; 00870 case 5: 00871 if(indexar[0]==0) accum_mballfaces(indexar[0], indexar[3], indexar[2], indexar[1]); 00872 else accum_mballfaces(indexar[3], indexar[2], indexar[1], indexar[0]); 00873 00874 accum_mballfaces(indexar[4], indexar[3], indexar[0], 0); 00875 break; 00876 case 6: 00877 if(indexar[0]==0) { 00878 accum_mballfaces(indexar[0], indexar[3], indexar[2], indexar[1]); 00879 accum_mballfaces(indexar[0], indexar[5], indexar[4], indexar[3]); 00880 } 00881 else { 00882 accum_mballfaces(indexar[3], indexar[2], indexar[1], indexar[0]); 00883 accum_mballfaces(indexar[5], indexar[4], indexar[3], indexar[0]); 00884 } 00885 break; 00886 case 7: 00887 if(indexar[0]==0) { 00888 accum_mballfaces(indexar[0], indexar[3], indexar[2], indexar[1]); 00889 accum_mballfaces(indexar[0], indexar[5], indexar[4], indexar[3]); 00890 } 00891 else { 00892 accum_mballfaces(indexar[3], indexar[2], indexar[1], indexar[0]); 00893 accum_mballfaces(indexar[5], indexar[4], indexar[3], indexar[0]); 00894 } 00895 00896 accum_mballfaces(indexar[6], indexar[5], indexar[0], 0); 00897 00898 break; 00899 } 00900 } 00901 } 00902 } 00903 00904 00905 /* testface: given cube at lattice (i, j, k), and four corners of face, 00906 * if surface crosses face, compute other four corners of adjacent cube 00907 * and add new cube to cube stack */ 00908 00909 void testface(int i, int j, int k, CUBE* old, int bit, int c1, int c2, int c3, int c4, PROCESS *p) 00910 { 00911 CUBE newc; 00912 CUBES *oldcubes = p->cubes; 00913 CORNER *corn1, *corn2, *corn3, *corn4; 00914 int n, pos; 00915 00916 corn1= old->corners[c1]; 00917 corn2= old->corners[c2]; 00918 corn3= old->corners[c3]; 00919 corn4= old->corners[c4]; 00920 00921 pos = corn1->value > 0.0f ? 1 : 0; 00922 00923 /* test if no surface crossing */ 00924 if( (corn2->value > 0) == pos && (corn3->value > 0) == pos && (corn4->value > 0) == pos) return; 00925 /* test if cube out of bounds */ 00926 /*if ( abs(i) > p->bounds || abs(j) > p->bounds || abs(k) > p->bounds) return;*/ 00927 /* test if already visited (always as last) */ 00928 if (setcenter(p->centers, i, j, k)) return; 00929 00930 00931 /* create new cube and add cube to top of stack: */ 00932 p->cubes = (CUBES *) new_pgn_element(sizeof(CUBES)); 00933 p->cubes->next = oldcubes; 00934 00935 newc.i = i; 00936 newc.j = j; 00937 newc.k = k; 00938 for (n = 0; n < 8; n++) newc.corners[n] = NULL; 00939 00940 newc.corners[FLIP(c1, bit)] = corn1; 00941 newc.corners[FLIP(c2, bit)] = corn2; 00942 newc.corners[FLIP(c3, bit)] = corn3; 00943 newc.corners[FLIP(c4, bit)] = corn4; 00944 00945 if(newc.corners[0]==NULL) newc.corners[0] = setcorner(p, i, j, k); 00946 if(newc.corners[1]==NULL) newc.corners[1] = setcorner(p, i, j, k+1); 00947 if(newc.corners[2]==NULL) newc.corners[2] = setcorner(p, i, j+1, k); 00948 if(newc.corners[3]==NULL) newc.corners[3] = setcorner(p, i, j+1, k+1); 00949 if(newc.corners[4]==NULL) newc.corners[4] = setcorner(p, i+1, j, k); 00950 if(newc.corners[5]==NULL) newc.corners[5] = setcorner(p, i+1, j, k+1); 00951 if(newc.corners[6]==NULL) newc.corners[6] = setcorner(p, i+1, j+1, k); 00952 if(newc.corners[7]==NULL) newc.corners[7] = setcorner(p, i+1, j+1, k+1); 00953 00954 p->cubes->cube= newc; 00955 } 00956 00957 /* setcorner: return corner with the given lattice location 00958 set (and cache) its function value */ 00959 00960 CORNER *setcorner (PROCESS* p, int i, int j, int k) 00961 { 00962 /* for speed, do corner value caching here */ 00963 CORNER *c; 00964 int index; 00965 00966 /* does corner exist? */ 00967 index = HASH(i, j, k); 00968 c = p->corners[index]; 00969 00970 for (; c != NULL; c = c->next) { 00971 if (c->i == i && c->j == j && c->k == k) { 00972 return c; 00973 } 00974 } 00975 00976 c = (CORNER *) new_pgn_element(sizeof(CORNER)); 00977 00978 c->i = i; 00979 c->x = ((float)i-0.5f)*p->size; 00980 c->j = j; 00981 c->y = ((float)j-0.5f)*p->size; 00982 c->k = k; 00983 c->z = ((float)k-0.5f)*p->size; 00984 c->value = p->function(c->x, c->y, c->z); 00985 00986 c->next = p->corners[index]; 00987 p->corners[index] = c; 00988 00989 return c; 00990 } 00991 00992 00993 /* nextcwedge: return next clockwise edge from given edge around given face */ 00994 00995 int nextcwedge (int edge, int face) 00996 { 00997 switch (edge) { 00998 case LB: 00999 return (face == L)? LF : BN; 01000 case LT: 01001 return (face == L)? LN : TF; 01002 case LN: 01003 return (face == L)? LB : TN; 01004 case LF: 01005 return (face == L)? LT : BF; 01006 case RB: 01007 return (face == R)? RN : BF; 01008 case RT: 01009 return (face == R)? RF : TN; 01010 case RN: 01011 return (face == R)? RT : BN; 01012 case RF: 01013 return (face == R)? RB : TF; 01014 case BN: 01015 return (face == B)? RB : LN; 01016 case BF: 01017 return (face == B)? LB : RF; 01018 case TN: 01019 return (face == T)? LT : RN; 01020 case TF: 01021 return (face == T)? RT : LF; 01022 } 01023 return 0; 01024 } 01025 01026 01027 /* otherface: return face adjoining edge that is not the given face */ 01028 01029 int otherface (int edge, int face) 01030 { 01031 int other = leftface[edge]; 01032 return face == other? rightface[edge] : other; 01033 } 01034 01035 01036 /* makecubetable: create the 256 entry table for cubical polygonization */ 01037 01038 void makecubetable (void) 01039 { 01040 static int isdone= 0; 01041 int i, e, c, done[12], pos[8]; 01042 01043 if(isdone) return; 01044 isdone= 1; 01045 01046 for (i = 0; i < 256; i++) { 01047 for (e = 0; e < 12; e++) done[e] = 0; 01048 for (c = 0; c < 8; c++) pos[c] = MB_BIT(i, c); 01049 for (e = 0; e < 12; e++) 01050 if (!done[e] && (pos[corner1[e]] != pos[corner2[e]])) { 01051 INTLIST *ints = NULL; 01052 INTLISTS *lists = (INTLISTS *) MEM_callocN(sizeof(INTLISTS), "mball_intlist"); 01053 int start = e, edge = e; 01054 01055 /* get face that is to right of edge from pos to neg corner: */ 01056 int face = pos[corner1[e]]? rightface[e] : leftface[e]; 01057 01058 while (1) { 01059 edge = nextcwedge(edge, face); 01060 done[edge] = 1; 01061 if (pos[corner1[edge]] != pos[corner2[edge]]) { 01062 INTLIST *tmp = ints; 01063 01064 ints = (INTLIST *) MEM_callocN(sizeof(INTLIST), "mball_intlist"); 01065 ints->i = edge; 01066 ints->next = tmp; /* add edge to head of list */ 01067 01068 if (edge == start) break; 01069 face = otherface(edge, face); 01070 } 01071 } 01072 lists->list = ints; /* add ints to head of table entry */ 01073 lists->next = cubetable[i]; 01074 cubetable[i] = lists; 01075 } 01076 } 01077 } 01078 01079 void BKE_freecubetable(void) 01080 { 01081 int i; 01082 INTLISTS *lists, *nlists; 01083 INTLIST *ints, *nints; 01084 01085 for (i = 0; i < 256; i++) { 01086 lists= cubetable[i]; 01087 while(lists) { 01088 nlists= lists->next; 01089 01090 ints= lists->list; 01091 while(ints) { 01092 nints= ints->next; 01093 MEM_freeN(ints); 01094 ints= nints; 01095 } 01096 01097 MEM_freeN(lists); 01098 lists= nlists; 01099 } 01100 cubetable[i]= NULL; 01101 } 01102 } 01103 01104 /**** Storage ****/ 01105 01106 /* setcenter: set (i,j,k) entry of table[] 01107 * return 1 if already set; otherwise, set and return 0 */ 01108 01109 int setcenter(CENTERLIST *table[], int i, int j, int k) 01110 { 01111 int index; 01112 CENTERLIST *newc, *l, *q; 01113 01114 index= HASH(i, j, k); 01115 q= table[index]; 01116 01117 for (l = q; l != NULL; l = l->next) { 01118 if (l->i == i && l->j == j && l->k == k) return 1; 01119 } 01120 01121 newc = (CENTERLIST *) new_pgn_element(sizeof(CENTERLIST)); 01122 newc->i = i; 01123 newc->j = j; 01124 newc->k = k; 01125 newc->next = q; 01126 table[index] = newc; 01127 01128 return 0; 01129 } 01130 01131 01132 /* setedge: set vertex id for edge */ 01133 01134 void setedge (EDGELIST *table[], 01135 int i1, int j1, 01136 int k1, int i2, 01137 int j2, int k2, 01138 int vid) 01139 { 01140 unsigned int index; 01141 EDGELIST *newe; 01142 01143 if (i1>i2 || (i1==i2 && (j1>j2 || (j1==j2 && k1>k2)))) { 01144 int t=i1; 01145 i1=i2; 01146 i2=t; 01147 t=j1; 01148 j1=j2; 01149 j2=t; 01150 t=k1; 01151 k1=k2; 01152 k2=t; 01153 } 01154 index = HASH(i1, j1, k1) + HASH(i2, j2, k2); 01155 newe = (EDGELIST *) new_pgn_element(sizeof(EDGELIST)); 01156 newe->i1 = i1; 01157 newe->j1 = j1; 01158 newe->k1 = k1; 01159 newe->i2 = i2; 01160 newe->j2 = j2; 01161 newe->k2 = k2; 01162 newe->vid = vid; 01163 newe->next = table[index]; 01164 table[index] = newe; 01165 } 01166 01167 01168 /* getedge: return vertex id for edge; return -1 if not set */ 01169 01170 int getedge (EDGELIST *table[], 01171 int i1, int j1, int k1, 01172 int i2, int j2, int k2) 01173 { 01174 EDGELIST *q; 01175 01176 if (i1>i2 || (i1==i2 && (j1>j2 || (j1==j2 && k1>k2)))) { 01177 int t=i1; 01178 i1=i2; 01179 i2=t; 01180 t=j1; 01181 j1=j2; 01182 j2=t; 01183 t=k1; 01184 k1=k2; 01185 k2=t; 01186 } 01187 q = table[HASH(i1, j1, k1)+HASH(i2, j2, k2)]; 01188 for (; q != NULL; q = q->next) 01189 if (q->i1 == i1 && q->j1 == j1 && q->k1 == k1 && 01190 q->i2 == i2 && q->j2 == j2 && q->k2 == k2) 01191 return q->vid; 01192 return -1; 01193 } 01194 01195 01196 /**** Vertices ****/ 01197 01198 #undef R 01199 01200 01201 01202 /* vertid: return index for vertex on edge: 01203 * c1->value and c2->value are presumed of different sign 01204 * return saved index if any; else compute vertex and save */ 01205 01206 /* addtovertices: add v to sequence of vertices */ 01207 01208 void addtovertices (VERTICES *vertices, VERTEX v) 01209 { 01210 if (vertices->count == vertices->max) { 01211 int i; 01212 VERTEX *newv; 01213 vertices->max = vertices->count == 0 ? 10 : 2*vertices->count; 01214 newv = (VERTEX *) MEM_callocN(vertices->max * sizeof(VERTEX), "addtovertices"); 01215 01216 for (i = 0; i < vertices->count; i++) newv[i] = vertices->ptr[i]; 01217 01218 if (vertices->ptr != NULL) MEM_freeN(vertices->ptr); 01219 vertices->ptr = newv; 01220 } 01221 vertices->ptr[vertices->count++] = v; 01222 } 01223 01224 /* vnormal: compute unit length surface normal at point */ 01225 01226 void vnormal (MB_POINT *point, PROCESS *p, MB_POINT *v) 01227 { 01228 float delta= 0.2f*p->delta; 01229 float f = p->function(point->x, point->y, point->z); 01230 01231 v->x = p->function(point->x+delta, point->y, point->z)-f; 01232 v->y = p->function(point->x, point->y+delta, point->z)-f; 01233 v->z = p->function(point->x, point->y, point->z+delta)-f; 01234 f = sqrtf(v->x*v->x + v->y*v->y + v->z*v->z); 01235 01236 if (f != 0.0f) { 01237 v->x /= f; 01238 v->y /= f; 01239 v->z /= f; 01240 } 01241 01242 if(FALSE) { 01243 MB_POINT temp; 01244 01245 delta *= 2.0f; 01246 01247 f = p->function(point->x, point->y, point->z); 01248 01249 temp.x = p->function(point->x+delta, point->y, point->z)-f; 01250 temp.y = p->function(point->x, point->y+delta, point->z)-f; 01251 temp.z = p->function(point->x, point->y, point->z+delta)-f; 01252 f = sqrtf(temp.x*temp.x + temp.y*temp.y + temp.z*temp.z); 01253 01254 if (f != 0.0f) { 01255 temp.x /= f; 01256 temp.y /= f; 01257 temp.z /= f; 01258 01259 v->x+= temp.x; 01260 v->y+= temp.y; 01261 v->z+= temp.z; 01262 01263 f = sqrtf(v->x*v->x + v->y*v->y + v->z*v->z); 01264 01265 if (f != 0.0f) { 01266 v->x /= f; 01267 v->y /= f; 01268 v->z /= f; 01269 } 01270 } 01271 } 01272 01273 } 01274 01275 01276 int vertid (CORNER *c1, CORNER *c2, PROCESS *p, MetaBall *mb) 01277 { 01278 VERTEX v; 01279 MB_POINT a, b; 01280 int vid = getedge(p->edges, c1->i, c1->j, c1->k, c2->i, c2->j, c2->k); 01281 01282 if (vid != -1) return vid; /* previously computed */ 01283 a.x = c1->x; 01284 a.y = c1->y; 01285 a.z = c1->z; 01286 b.x = c2->x; 01287 b.y = c2->y; 01288 b.z = c2->z; 01289 01290 converge(&a, &b, c1->value, c2->value, p->function, &v.position, mb, 1); /* position */ 01291 vnormal(&v.position, p, &v.normal); 01292 01293 addtovertices(&p->vertices, v); /* save vertex */ 01294 vid = p->vertices.count-1; 01295 setedge(p->edges, c1->i, c1->j, c1->k, c2->i, c2->j, c2->k, vid); 01296 01297 return vid; 01298 } 01299 01300 01301 01302 01303 /* converge: from two points of differing sign, converge to zero crossing */ 01304 /* watch it: p1 and p2 are used to calculate */ 01305 void converge (MB_POINT *p1, MB_POINT *p2, float v1, float v2, 01306 float (*function)(float, float, float), MB_POINT *p, MetaBall *mb, int f) 01307 { 01308 int i = 0; 01309 MB_POINT pos, neg; 01310 float positive = 0.0f, negative = 0.0f; 01311 float dx = 0.0f ,dy = 0.0f ,dz = 0.0f; 01312 01313 if (v1 < 0) { 01314 pos= *p2; 01315 neg= *p1; 01316 positive = v2; 01317 negative = v1; 01318 } 01319 else { 01320 pos= *p1; 01321 neg= *p2; 01322 positive = v1; 01323 negative = v2; 01324 } 01325 01326 dx = pos.x - neg.x; 01327 dy = pos.y - neg.y; 01328 dz = pos.z - neg.z; 01329 01330 /* Approximation by linear interpolation is faster then binary subdivision, 01331 * but it results sometimes (mb->thresh < 0.2) into the strange results */ 01332 if((mb->thresh > 0.2f) && (f==1)){ 01333 if((dy == 0.0f) && (dz == 0.0f)){ 01334 p->x = neg.x - negative*dx/(positive-negative); 01335 p->y = neg.y; 01336 p->z = neg.z; 01337 return; 01338 } 01339 if((dx == 0.0f) && (dz == 0.0f)){ 01340 p->x = neg.x; 01341 p->y = neg.y - negative*dy/(positive-negative); 01342 p->z = neg.z; 01343 return; 01344 } 01345 if((dx == 0.0f) && (dy == 0.0f)){ 01346 p->x = neg.x; 01347 p->y = neg.y; 01348 p->z = neg.z - negative*dz/(positive-negative); 01349 return; 01350 } 01351 } 01352 01353 if((dy == 0.0f) && (dz == 0.0f)){ 01354 p->y = neg.y; 01355 p->z = neg.z; 01356 while (1) { 01357 if (i++ == RES) return; 01358 p->x = 0.5f*(pos.x + neg.x); 01359 if ((function(p->x,p->y,p->z)) > 0.0f) pos.x = p->x; else neg.x = p->x; 01360 } 01361 } 01362 01363 if((dx == 0.0f) && (dz == 0.0f)){ 01364 p->x = neg.x; 01365 p->z = neg.z; 01366 while (1) { 01367 if (i++ == RES) return; 01368 p->y = 0.5f*(pos.y + neg.y); 01369 if ((function(p->x,p->y,p->z)) > 0.0f) pos.y = p->y; else neg.y = p->y; 01370 } 01371 } 01372 01373 if((dx == 0.0f) && (dy == 0.0f)){ 01374 p->x = neg.x; 01375 p->y = neg.y; 01376 while (1) { 01377 if (i++ == RES) return; 01378 p->z = 0.5f*(pos.z + neg.z); 01379 if ((function(p->x,p->y,p->z)) > 0.0f) pos.z = p->z; else neg.z = p->z; 01380 } 01381 } 01382 01383 /* This is necessary to find start point */ 01384 while (1) { 01385 p->x = 0.5f*(pos.x + neg.x); 01386 p->y = 0.5f*(pos.y + neg.y); 01387 p->z = 0.5f*(pos.z + neg.z); 01388 01389 if (i++ == RES) return; 01390 01391 if ((function(p->x, p->y, p->z)) > 0.0f){ 01392 pos.x = p->x; 01393 pos.y = p->y; 01394 pos.z = p->z; 01395 } 01396 else{ 01397 neg.x = p->x; 01398 neg.y = p->y; 01399 neg.z = p->z; 01400 } 01401 } 01402 } 01403 01404 /* ************************************** */ 01405 void add_cube(PROCESS *mbproc, int i, int j, int k, int count) 01406 { 01407 CUBES *ncube; 01408 int n; 01409 int a, b, c; 01410 01411 /* hmmm, not only one, but eight cube will be added on the stack 01412 * ... */ 01413 for(a=i-1; a<i+count; a++) 01414 for(b=j-1; b<j+count; b++) 01415 for(c=k-1; c<k+count; c++) { 01416 /* test if cube has been found before */ 01417 if( setcenter(mbproc->centers, a, b, c)==0 ) { 01418 /* push cube on stack: */ 01419 ncube= (CUBES *) new_pgn_element(sizeof(CUBES)); 01420 ncube->next= mbproc->cubes; 01421 mbproc->cubes= ncube; 01422 01423 ncube->cube.i= a; 01424 ncube->cube.j= b; 01425 ncube->cube.k= c; 01426 01427 /* set corners of initial cube: */ 01428 for (n = 0; n < 8; n++) 01429 ncube->cube.corners[n] = setcorner(mbproc, a+MB_BIT(n,2), b+MB_BIT(n,1), c+MB_BIT(n,0)); 01430 } 01431 } 01432 } 01433 01434 01435 void find_first_points(PROCESS *mbproc, MetaBall *mb, int a) 01436 { 01437 MB_POINT IN, in, OUT, out; /*point;*/ 01438 MetaElem *ml; 01439 int i, j, k, c_i, c_j, c_k; 01440 int index[3]={1,0,-1}; 01441 float f =0.0f; 01442 float in_v /*, out_v*/; 01443 MB_POINT workp; 01444 float tmp_v, workp_v, max_len, len, dx, dy, dz, nx, ny, nz, MAXN; 01445 01446 ml = mainb[a]; 01447 01448 f = 1-(mb->thresh/ml->s); 01449 01450 /* Skip, when Stiffness of MetaElement is too small ... MetaElement can't be 01451 * visible alone ... but still can influence others MetaElements :-) */ 01452 if(f > 0.0f) { 01453 OUT.x = IN.x = in.x= 0.0; 01454 OUT.y = IN.y = in.y= 0.0; 01455 OUT.z = IN.z = in.z= 0.0; 01456 01457 calc_mballco(ml, (float *)&in); 01458 in_v = mbproc->function(in.x, in.y, in.z); 01459 01460 for(i=0;i<3;i++){ 01461 switch (ml->type) { 01462 case MB_BALL: 01463 OUT.x = out.x= IN.x + index[i]*ml->rad; 01464 break; 01465 case MB_TUBE: 01466 case MB_PLANE: 01467 case MB_ELIPSOID: 01468 case MB_CUBE: 01469 OUT.x = out.x= IN.x + index[i]*(ml->expx + ml->rad); 01470 break; 01471 } 01472 01473 for(j=0;j<3;j++) { 01474 switch (ml->type) { 01475 case MB_BALL: 01476 OUT.y = out.y= IN.y + index[j]*ml->rad; 01477 break; 01478 case MB_TUBE: 01479 case MB_PLANE: 01480 case MB_ELIPSOID: 01481 case MB_CUBE: 01482 OUT.y = out.y= IN.y + index[j]*(ml->expy + ml->rad); 01483 break; 01484 } 01485 01486 for(k=0;k<3;k++) { 01487 out.x = OUT.x; 01488 out.y = OUT.y; 01489 switch (ml->type) { 01490 case MB_BALL: 01491 case MB_TUBE: 01492 case MB_PLANE: 01493 out.z= IN.z + index[k]*ml->rad; 01494 break; 01495 case MB_ELIPSOID: 01496 case MB_CUBE: 01497 out.z= IN.z + index[k]*(ml->expz + ml->rad); 01498 break; 01499 } 01500 01501 calc_mballco(ml, (float *)&out); 01502 01503 /*out_v = mbproc->function(out.x, out.y, out.z);*/ /*UNUSED*/ 01504 01505 /* find "first points" on Implicit Surface of MetaElemnt ml */ 01506 workp.x = in.x; 01507 workp.y = in.y; 01508 workp.z = in.z; 01509 workp_v = in_v; 01510 max_len = sqrtf((out.x-in.x)*(out.x-in.x) + (out.y-in.y)*(out.y-in.y) + (out.z-in.z)*(out.z-in.z)); 01511 01512 nx = abs((out.x - in.x)/mbproc->size); 01513 ny = abs((out.y - in.y)/mbproc->size); 01514 nz = abs((out.z - in.z)/mbproc->size); 01515 01516 MAXN = MAX3(nx,ny,nz); 01517 if(MAXN!=0.0f) { 01518 dx = (out.x - in.x)/MAXN; 01519 dy = (out.y - in.y)/MAXN; 01520 dz = (out.z - in.z)/MAXN; 01521 01522 len = 0.0; 01523 while(len<=max_len) { 01524 workp.x += dx; 01525 workp.y += dy; 01526 workp.z += dz; 01527 /* compute value of implicite function */ 01528 tmp_v = mbproc->function(workp.x, workp.y, workp.z); 01529 /* add cube to the stack, when value of implicite function crosses zero value */ 01530 if((tmp_v<0.0f && workp_v>=0.0f)||(tmp_v>0.0f && workp_v<=0.0f)) { 01531 01532 /* indexes of CUBE, which includes "first point" */ 01533 c_i= (int)floor(workp.x/mbproc->size); 01534 c_j= (int)floor(workp.y/mbproc->size); 01535 c_k= (int)floor(workp.z/mbproc->size); 01536 01537 /* add CUBE (with indexes c_i, c_j, c_k) to the stack, 01538 * this cube includes found point of Implicit Surface */ 01539 if (ml->flag & MB_NEGATIVE) 01540 add_cube(mbproc, c_i, c_j, c_k, 2); 01541 else 01542 add_cube(mbproc, c_i, c_j, c_k, 1); 01543 } 01544 len = sqrtf((workp.x-in.x)*(workp.x-in.x) + (workp.y-in.y)*(workp.y-in.y) + (workp.z-in.z)*(workp.z-in.z)); 01545 workp_v = tmp_v; 01546 01547 } 01548 } 01549 } 01550 } 01551 } 01552 } 01553 } 01554 01555 void polygonize(PROCESS *mbproc, MetaBall *mb) 01556 { 01557 CUBE c; 01558 int a; 01559 01560 mbproc->vertices.count = mbproc->vertices.max = 0; 01561 mbproc->vertices.ptr = NULL; 01562 01563 /* allocate hash tables and build cube polygon table: */ 01564 mbproc->centers = MEM_callocN(HASHSIZE * sizeof(CENTERLIST *), "mbproc->centers"); 01565 mbproc->corners = MEM_callocN(HASHSIZE * sizeof(CORNER *), "mbproc->corners"); 01566 mbproc->edges = MEM_callocN(2*HASHSIZE * sizeof(EDGELIST *), "mbproc->edges"); 01567 makecubetable(); 01568 01569 for(a=0; a<totelem; a++) { 01570 01571 /* try to find 8 points on the surface for each MetaElem */ 01572 find_first_points(mbproc, mb, a); 01573 } 01574 01575 /* polygonize all MetaElems of current MetaBall */ 01576 while (mbproc->cubes != NULL) { /* process active cubes till none left */ 01577 c = mbproc->cubes->cube; 01578 01579 /* polygonize the cube directly: */ 01580 docube(&c, mbproc, mb); 01581 01582 /* pop current cube from stack */ 01583 mbproc->cubes = mbproc->cubes->next; 01584 01585 /* test six face directions, maybe add to stack: */ 01586 testface(c.i-1, c.j, c.k, &c, 2, LBN, LBF, LTN, LTF, mbproc); 01587 testface(c.i+1, c.j, c.k, &c, 2, RBN, RBF, RTN, RTF, mbproc); 01588 testface(c.i, c.j-1, c.k, &c, 1, LBN, LBF, RBN, RBF, mbproc); 01589 testface(c.i, c.j+1, c.k, &c, 1, LTN, LTF, RTN, RTF, mbproc); 01590 testface(c.i, c.j, c.k-1, &c, 0, LBN, LTN, RBN, RTN, mbproc); 01591 testface(c.i, c.j, c.k+1, &c, 0, LBF, LTF, RBF, RTF, mbproc); 01592 } 01593 } 01594 01595 float init_meta(Scene *scene, Object *ob) /* return totsize */ 01596 { 01597 Scene *sce_iter= scene; 01598 Base *base; 01599 Object *bob; 01600 MetaBall *mb; 01601 MetaElem *ml; 01602 float size, totsize, obinv[4][4], obmat[4][4], vec[3]; 01603 //float max=0.0; 01604 int a, obnr, zero_size=0; 01605 char obname[32]; 01606 01607 copy_m4_m4(obmat, ob->obmat); /* to cope with duplicators from next_object */ 01608 invert_m4_m4(obinv, ob->obmat); 01609 a= 0; 01610 01611 BLI_split_name_num(obname, &obnr, ob->id.name+2, '.'); 01612 01613 /* make main array */ 01614 next_object(&sce_iter, 0, NULL, NULL); 01615 while(next_object(&sce_iter, 1, &base, &bob)) { 01616 01617 if(bob->type==OB_MBALL) { 01618 zero_size= 0; 01619 ml= NULL; 01620 01621 if(bob==ob && (base->flag & OB_FROMDUPLI)==0) { 01622 mb= ob->data; 01623 01624 if(mb->editelems) ml= mb->editelems->first; 01625 else ml= mb->elems.first; 01626 } 01627 else { 01628 char name[32]; 01629 int nr; 01630 01631 BLI_split_name_num(name, &nr, bob->id.name+2, '.'); 01632 if( strcmp(obname, name)==0 ) { 01633 mb= bob->data; 01634 01635 if(mb->editelems) ml= mb->editelems->first; 01636 else ml= mb->elems.first; 01637 } 01638 } 01639 01640 /* when metaball object has zero scale, then MetaElem to this MetaBall 01641 * will not be put to mainb array */ 01642 if(bob->size[0]==0.0f || bob->size[1]==0.0f || bob->size[2]==0.0f) { 01643 zero_size= 1; 01644 } 01645 else if(bob->parent) { 01646 struct Object *pob=bob->parent; 01647 while(pob) { 01648 if(pob->size[0]==0.0f || pob->size[1]==0.0f || pob->size[2]==0.0f) { 01649 zero_size= 1; 01650 break; 01651 } 01652 pob= pob->parent; 01653 } 01654 } 01655 01656 if (zero_size) { 01657 unsigned int ml_count=0; 01658 while(ml) { 01659 ml_count++; 01660 ml= ml->next; 01661 } 01662 totelem -= ml_count; 01663 } 01664 else { 01665 while(ml) { 01666 if(!(ml->flag & MB_HIDE)) { 01667 int i; 01668 float temp1[4][4], temp2[4][4], temp3[4][4]; 01669 float (*mat)[4] = NULL, (*imat)[4] = NULL; 01670 float max_x, max_y, max_z, min_x, min_y, min_z; 01671 01672 max_x = max_y = max_z = -3.4e38; 01673 min_x = min_y = min_z = 3.4e38; 01674 01675 /* too big stiffness seems only ugly due to linear interpolation 01676 * no need to have possibility for too big stiffness */ 01677 if(ml->s > 10.0f) ml->s = 10.0f; 01678 01679 /* Rotation of MetaElem is stored in quat */ 01680 quat_to_mat4( temp3,ml->quat); 01681 01682 /* Translation of MetaElem */ 01683 unit_m4(temp2); 01684 temp2[3][0]= ml->x; 01685 temp2[3][1]= ml->y; 01686 temp2[3][2]= ml->z; 01687 01688 mul_m4_m4m4(temp1, temp3, temp2); 01689 01690 /* make a copy because of duplicates */ 01691 mainb[a]= new_pgn_element(sizeof(MetaElem)); 01692 *(mainb[a])= *ml; 01693 mainb[a]->bb = new_pgn_element(sizeof(BoundBox)); 01694 01695 mat= new_pgn_element(4*4*sizeof(float)); 01696 imat= new_pgn_element(4*4*sizeof(float)); 01697 01698 /* mat is the matrix to transform from mball into the basis-mball */ 01699 invert_m4_m4(obinv, obmat); 01700 mul_m4_m4m4(temp2, bob->obmat, obinv); 01701 /* MetaBall transformation */ 01702 mul_m4_m4m4(mat, temp1, temp2); 01703 01704 invert_m4_m4(imat,mat); 01705 01706 mainb[a]->rad2= ml->rad*ml->rad; 01707 01708 mainb[a]->mat= (float*) mat; 01709 mainb[a]->imat= (float*) imat; 01710 01711 /* untransformed Bounding Box of MetaElem */ 01712 /* 0 */ 01713 mainb[a]->bb->vec[0][0]= -ml->expx; 01714 mainb[a]->bb->vec[0][1]= -ml->expy; 01715 mainb[a]->bb->vec[0][2]= -ml->expz; 01716 /* 1 */ 01717 mainb[a]->bb->vec[1][0]= ml->expx; 01718 mainb[a]->bb->vec[1][1]= -ml->expy; 01719 mainb[a]->bb->vec[1][2]= -ml->expz; 01720 /* 2 */ 01721 mainb[a]->bb->vec[2][0]= ml->expx; 01722 mainb[a]->bb->vec[2][1]= ml->expy; 01723 mainb[a]->bb->vec[2][2]= -ml->expz; 01724 /* 3 */ 01725 mainb[a]->bb->vec[3][0]= -ml->expx; 01726 mainb[a]->bb->vec[3][1]= ml->expy; 01727 mainb[a]->bb->vec[3][2]= -ml->expz; 01728 /* 4 */ 01729 mainb[a]->bb->vec[4][0]= -ml->expx; 01730 mainb[a]->bb->vec[4][1]= -ml->expy; 01731 mainb[a]->bb->vec[4][2]= ml->expz; 01732 /* 5 */ 01733 mainb[a]->bb->vec[5][0]= ml->expx; 01734 mainb[a]->bb->vec[5][1]= -ml->expy; 01735 mainb[a]->bb->vec[5][2]= ml->expz; 01736 /* 6 */ 01737 mainb[a]->bb->vec[6][0]= ml->expx; 01738 mainb[a]->bb->vec[6][1]= ml->expy; 01739 mainb[a]->bb->vec[6][2]= ml->expz; 01740 /* 7 */ 01741 mainb[a]->bb->vec[7][0]= -ml->expx; 01742 mainb[a]->bb->vec[7][1]= ml->expy; 01743 mainb[a]->bb->vec[7][2]= ml->expz; 01744 01745 /* transformation of Metalem bb */ 01746 for(i=0; i<8; i++) 01747 mul_m4_v3((float ( * )[4])mat, mainb[a]->bb->vec[i]); 01748 01749 /* find max and min of transformed bb */ 01750 for(i=0; i<8; i++){ 01751 /* find maximums */ 01752 if(mainb[a]->bb->vec[i][0] > max_x) max_x = mainb[a]->bb->vec[i][0]; 01753 if(mainb[a]->bb->vec[i][1] > max_y) max_y = mainb[a]->bb->vec[i][1]; 01754 if(mainb[a]->bb->vec[i][2] > max_z) max_z = mainb[a]->bb->vec[i][2]; 01755 /* find minimums */ 01756 if(mainb[a]->bb->vec[i][0] < min_x) min_x = mainb[a]->bb->vec[i][0]; 01757 if(mainb[a]->bb->vec[i][1] < min_y) min_y = mainb[a]->bb->vec[i][1]; 01758 if(mainb[a]->bb->vec[i][2] < min_z) min_z = mainb[a]->bb->vec[i][2]; 01759 } 01760 01761 /* create "new" bb, only point 0 and 6, which are 01762 * neccesary for octal tree filling */ 01763 mainb[a]->bb->vec[0][0] = min_x - ml->rad; 01764 mainb[a]->bb->vec[0][1] = min_y - ml->rad; 01765 mainb[a]->bb->vec[0][2] = min_z - ml->rad; 01766 01767 mainb[a]->bb->vec[6][0] = max_x + ml->rad; 01768 mainb[a]->bb->vec[6][1] = max_y + ml->rad; 01769 mainb[a]->bb->vec[6][2] = max_z + ml->rad; 01770 01771 a++; 01772 } 01773 ml= ml->next; 01774 } 01775 } 01776 } 01777 } 01778 01779 01780 /* totsize (= 'manhattan' radius) */ 01781 totsize= 0.0; 01782 for(a=0; a<totelem; a++) { 01783 01784 vec[0]= mainb[a]->x + mainb[a]->rad + mainb[a]->expx; 01785 vec[1]= mainb[a]->y + mainb[a]->rad + mainb[a]->expy; 01786 vec[2]= mainb[a]->z + mainb[a]->rad + mainb[a]->expz; 01787 01788 calc_mballco(mainb[a], vec); 01789 01790 size= (float)fabs( vec[0] ); 01791 if( size > totsize ) totsize= size; 01792 size= (float)fabs( vec[1] ); 01793 if( size > totsize ) totsize= size; 01794 size= (float)fabs( vec[2] ); 01795 if( size > totsize ) totsize= size; 01796 01797 vec[0]= mainb[a]->x - mainb[a]->rad; 01798 vec[1]= mainb[a]->y - mainb[a]->rad; 01799 vec[2]= mainb[a]->z - mainb[a]->rad; 01800 01801 calc_mballco(mainb[a], vec); 01802 01803 size= (float)fabs( vec[0] ); 01804 if( size > totsize ) totsize= size; 01805 size= (float)fabs( vec[1] ); 01806 if( size > totsize ) totsize= size; 01807 size= (float)fabs( vec[2] ); 01808 if( size > totsize ) totsize= size; 01809 } 01810 01811 for(a=0; a<totelem; a++) { 01812 thresh+= densfunc( mainb[a], 2.0f*totsize, 2.0f*totsize, 2.0f*totsize); 01813 } 01814 01815 return totsize; 01816 } 01817 01818 /* if MetaElem lies in node, then node includes MetaElem pointer (ml_p) 01819 * pointing at MetaElem (ml) 01820 */ 01821 void fill_metaball_octal_node(octal_node *node, MetaElem *ml, short i) 01822 { 01823 ml_pointer *ml_p; 01824 01825 ml_p= MEM_mallocN(sizeof(ml_pointer), "ml_pointer"); 01826 ml_p->ml= ml; 01827 BLI_addtail(&(node->nodes[i]->elems), ml_p); 01828 node->count++; 01829 01830 if(ml->flag & MB_NEGATIVE) { 01831 node->nodes[i]->neg++; 01832 } 01833 else{ 01834 node->nodes[i]->pos++; 01835 } 01836 } 01837 01838 /* Node is subdivided as is ilustrated on the following figure: 01839 * 01840 * +------+------+ 01841 * / / /| 01842 * +------+------+ | 01843 * / / /| + 01844 * +------+------+ |/| 01845 * | | | + | 01846 * | | |/| + 01847 * +------+------+ |/ 01848 * | | | + 01849 * | | |/ 01850 * +------+------+ 01851 * 01852 */ 01853 void subdivide_metaball_octal_node(octal_node *node, float size_x, float size_y, float size_z, short depth) 01854 { 01855 MetaElem *ml; 01856 ml_pointer *ml_p; 01857 float x,y,z; 01858 int a,i; 01859 01860 /* create new nodes */ 01861 for(a=0;a<8;a++){ 01862 node->nodes[a]= MEM_mallocN(sizeof(octal_node),"octal_node"); 01863 for(i=0;i<8;i++) 01864 node->nodes[a]->nodes[i]= NULL; 01865 node->nodes[a]->parent= node; 01866 node->nodes[a]->elems.first= NULL; 01867 node->nodes[a]->elems.last= NULL; 01868 node->nodes[a]->count= 0; 01869 node->nodes[a]->neg= 0; 01870 node->nodes[a]->pos= 0; 01871 } 01872 01873 size_x /= 2; 01874 size_y /= 2; 01875 size_z /= 2; 01876 01877 /* center of node */ 01878 node->x = x = node->x_min + size_x; 01879 node->y = y = node->y_min + size_y; 01880 node->z = z = node->z_min + size_z; 01881 01882 /* setting up of border points of new nodes */ 01883 node->nodes[0]->x_min = node->x_min; 01884 node->nodes[0]->y_min = node->y_min; 01885 node->nodes[0]->z_min = node->z_min; 01886 node->nodes[0]->x = node->nodes[0]->x_min + size_x/2; 01887 node->nodes[0]->y = node->nodes[0]->y_min + size_y/2; 01888 node->nodes[0]->z = node->nodes[0]->z_min + size_z/2; 01889 01890 node->nodes[1]->x_min = x; 01891 node->nodes[1]->y_min = node->y_min; 01892 node->nodes[1]->z_min = node->z_min; 01893 node->nodes[1]->x = node->nodes[1]->x_min + size_x/2; 01894 node->nodes[1]->y = node->nodes[1]->y_min + size_y/2; 01895 node->nodes[1]->z = node->nodes[1]->z_min + size_z/2; 01896 01897 node->nodes[2]->x_min = x; 01898 node->nodes[2]->y_min = y; 01899 node->nodes[2]->z_min = node->z_min; 01900 node->nodes[2]->x = node->nodes[2]->x_min + size_x/2; 01901 node->nodes[2]->y = node->nodes[2]->y_min + size_y/2; 01902 node->nodes[2]->z = node->nodes[2]->z_min + size_z/2; 01903 01904 node->nodes[3]->x_min = node->x_min; 01905 node->nodes[3]->y_min = y; 01906 node->nodes[3]->z_min = node->z_min; 01907 node->nodes[3]->x = node->nodes[3]->x_min + size_x/2; 01908 node->nodes[3]->y = node->nodes[3]->y_min + size_y/2; 01909 node->nodes[3]->z = node->nodes[3]->z_min + size_z/2; 01910 01911 node->nodes[4]->x_min = node->x_min; 01912 node->nodes[4]->y_min = node->y_min; 01913 node->nodes[4]->z_min = z; 01914 node->nodes[4]->x = node->nodes[4]->x_min + size_x/2; 01915 node->nodes[4]->y = node->nodes[4]->y_min + size_y/2; 01916 node->nodes[4]->z = node->nodes[4]->z_min + size_z/2; 01917 01918 node->nodes[5]->x_min = x; 01919 node->nodes[5]->y_min = node->y_min; 01920 node->nodes[5]->z_min = z; 01921 node->nodes[5]->x = node->nodes[5]->x_min + size_x/2; 01922 node->nodes[5]->y = node->nodes[5]->y_min + size_y/2; 01923 node->nodes[5]->z = node->nodes[5]->z_min + size_z/2; 01924 01925 node->nodes[6]->x_min = x; 01926 node->nodes[6]->y_min = y; 01927 node->nodes[6]->z_min = z; 01928 node->nodes[6]->x = node->nodes[6]->x_min + size_x/2; 01929 node->nodes[6]->y = node->nodes[6]->y_min + size_y/2; 01930 node->nodes[6]->z = node->nodes[6]->z_min + size_z/2; 01931 01932 node->nodes[7]->x_min = node->x_min; 01933 node->nodes[7]->y_min = y; 01934 node->nodes[7]->z_min = z; 01935 node->nodes[7]->x = node->nodes[7]->x_min + size_x/2; 01936 node->nodes[7]->y = node->nodes[7]->y_min + size_y/2; 01937 node->nodes[7]->z = node->nodes[7]->z_min + size_z/2; 01938 01939 ml_p= node->elems.first; 01940 01941 /* setting up references of MetaElems for new nodes */ 01942 while(ml_p){ 01943 ml= ml_p->ml; 01944 if(ml->bb->vec[0][2] < z){ 01945 if(ml->bb->vec[0][1] < y){ 01946 /* vec[0][0] lies in first octant */ 01947 if(ml->bb->vec[0][0] < x){ 01948 /* ml belongs to the (0)1st node */ 01949 fill_metaball_octal_node(node, ml, 0); 01950 01951 /* ml belongs to the (3)4th node */ 01952 if(ml->bb->vec[6][1] >= y){ 01953 fill_metaball_octal_node(node, ml, 3); 01954 01955 /* ml belongs to the (7)8th node */ 01956 if(ml->bb->vec[6][2] >= z){ 01957 fill_metaball_octal_node(node, ml, 7); 01958 } 01959 } 01960 01961 /* ml belongs to the (1)2nd node */ 01962 if(ml->bb->vec[6][0] >= x){ 01963 fill_metaball_octal_node(node, ml, 1); 01964 01965 /* ml belongs to the (5)6th node */ 01966 if(ml->bb->vec[6][2] >= z){ 01967 fill_metaball_octal_node(node, ml, 5); 01968 } 01969 } 01970 01971 /* ml belongs to the (2)3th node */ 01972 if((ml->bb->vec[6][0] >= x) && (ml->bb->vec[6][1] >= y)){ 01973 fill_metaball_octal_node(node, ml, 2); 01974 01975 /* ml belong to the (6)7th node */ 01976 if(ml->bb->vec[6][2] >= z){ 01977 fill_metaball_octal_node(node, ml, 6); 01978 } 01979 01980 } 01981 01982 /* ml belongs to the (4)5th node too */ 01983 if(ml->bb->vec[6][2] >= z){ 01984 fill_metaball_octal_node(node, ml, 4); 01985 } 01986 01987 01988 01989 } 01990 /* vec[0][0] is in the (1)second octant */ 01991 else{ 01992 /* ml belong to the (1)2nd node */ 01993 fill_metaball_octal_node(node, ml, 1); 01994 01995 /* ml belongs to the (2)3th node */ 01996 if(ml->bb->vec[6][1] >= y){ 01997 fill_metaball_octal_node(node, ml, 2); 01998 01999 /* ml belongs to the (6)7th node */ 02000 if(ml->bb->vec[6][2] >= z){ 02001 fill_metaball_octal_node(node, ml, 6); 02002 } 02003 02004 } 02005 02006 /* ml belongs to the (5)6th node */ 02007 if(ml->bb->vec[6][2] >= z){ 02008 fill_metaball_octal_node(node, ml, 5); 02009 } 02010 } 02011 } 02012 else{ 02013 /* vec[0][0] is in the (3)4th octant */ 02014 if(ml->bb->vec[0][0] < x){ 02015 /* ml belongs to the (3)4nd node */ 02016 fill_metaball_octal_node(node, ml, 3); 02017 02018 /* ml belongs to the (7)8th node */ 02019 if(ml->bb->vec[6][2] >= z){ 02020 fill_metaball_octal_node(node, ml, 7); 02021 } 02022 02023 02024 /* ml belongs to the (2)3th node */ 02025 if(ml->bb->vec[6][0] >= x){ 02026 fill_metaball_octal_node(node, ml, 2); 02027 02028 /* ml belongs to the (6)7th node */ 02029 if(ml->bb->vec[6][2] >= z){ 02030 fill_metaball_octal_node(node, ml, 6); 02031 } 02032 } 02033 } 02034 02035 } 02036 02037 /* vec[0][0] is in the (2)3th octant */ 02038 if((ml->bb->vec[0][0] >= x) && (ml->bb->vec[0][1] >= y)){ 02039 /* ml belongs to the (2)3th node */ 02040 fill_metaball_octal_node(node, ml, 2); 02041 02042 /* ml belongs to the (6)7th node */ 02043 if(ml->bb->vec[6][2] >= z){ 02044 fill_metaball_octal_node(node, ml, 6); 02045 } 02046 } 02047 } 02048 else{ 02049 if(ml->bb->vec[0][1] < y){ 02050 /* vec[0][0] lies in (4)5th octant */ 02051 if(ml->bb->vec[0][0] < x){ 02052 /* ml belongs to the (4)5th node */ 02053 fill_metaball_octal_node(node, ml, 4); 02054 02055 if(ml->bb->vec[6][0] >= x){ 02056 fill_metaball_octal_node(node, ml, 5); 02057 } 02058 02059 if(ml->bb->vec[6][1] >= y){ 02060 fill_metaball_octal_node(node, ml, 7); 02061 } 02062 02063 if((ml->bb->vec[6][0] >= x) && (ml->bb->vec[6][1] >= y)){ 02064 fill_metaball_octal_node(node, ml, 6); 02065 } 02066 } 02067 /* vec[0][0] lies in (5)6th octant */ 02068 else{ 02069 fill_metaball_octal_node(node, ml, 5); 02070 02071 if(ml->bb->vec[6][1] >= y){ 02072 fill_metaball_octal_node(node, ml, 6); 02073 } 02074 } 02075 } 02076 else{ 02077 /* vec[0][0] lies in (7)8th octant */ 02078 if(ml->bb->vec[0][0] < x){ 02079 fill_metaball_octal_node(node, ml, 7); 02080 02081 if(ml->bb->vec[6][0] >= x){ 02082 fill_metaball_octal_node(node, ml, 6); 02083 } 02084 } 02085 02086 } 02087 02088 /* vec[0][0] lies in (6)7th octant */ 02089 if((ml->bb->vec[0][0] >= x) && (ml->bb->vec[0][1] >= y)){ 02090 fill_metaball_octal_node(node, ml, 6); 02091 } 02092 } 02093 ml_p= ml_p->next; 02094 } 02095 02096 /* free references of MetaElems for curent node (it is not needed anymore) */ 02097 BLI_freelistN(&node->elems); 02098 02099 depth--; 02100 02101 if(depth>0){ 02102 for(a=0;a<8;a++){ 02103 if(node->nodes[a]->count > 0) /* if node is not empty, then it is subdivided */ 02104 subdivide_metaball_octal_node(node->nodes[a], size_x, size_y, size_z, depth); 02105 } 02106 } 02107 } 02108 02109 /* free all octal nodes recursively */ 02110 void free_metaball_octal_node(octal_node *node) 02111 { 02112 int a; 02113 for(a=0;a<8;a++){ 02114 if(node->nodes[a]!=NULL) free_metaball_octal_node(node->nodes[a]); 02115 } 02116 BLI_freelistN(&node->elems); 02117 MEM_freeN(node); 02118 } 02119 02120 /* If scene include more then one MetaElem, then octree is used */ 02121 void init_metaball_octal_tree(int depth) 02122 { 02123 struct octal_node *node; 02124 ml_pointer *ml_p; 02125 float size[3]; 02126 int a; 02127 02128 metaball_tree= MEM_mallocN(sizeof(octal_tree), "metaball_octal_tree"); 02129 metaball_tree->first= node= MEM_mallocN(sizeof(octal_node), "metaball_octal_node"); 02130 /* maximal depth of octree */ 02131 metaball_tree->depth= depth; 02132 02133 metaball_tree->neg= node->neg=0; 02134 metaball_tree->pos= node->pos=0; 02135 02136 node->elems.first= NULL; 02137 node->elems.last= NULL; 02138 node->count=0; 02139 02140 for(a=0;a<8;a++) 02141 node->nodes[a]=NULL; 02142 02143 node->x_min= node->y_min= node->z_min= FLT_MAX; 02144 node->x_max= node->y_max= node->z_max= -FLT_MAX; 02145 02146 /* size of octal tree scene */ 02147 for(a=0;a<totelem;a++) { 02148 if(mainb[a]->bb->vec[0][0] < node->x_min) node->x_min= mainb[a]->bb->vec[0][0]; 02149 if(mainb[a]->bb->vec[0][1] < node->y_min) node->y_min= mainb[a]->bb->vec[0][1]; 02150 if(mainb[a]->bb->vec[0][2] < node->z_min) node->z_min= mainb[a]->bb->vec[0][2]; 02151 02152 if(mainb[a]->bb->vec[6][0] > node->x_max) node->x_max= mainb[a]->bb->vec[6][0]; 02153 if(mainb[a]->bb->vec[6][1] > node->y_max) node->y_max= mainb[a]->bb->vec[6][1]; 02154 if(mainb[a]->bb->vec[6][2] > node->z_max) node->z_max= mainb[a]->bb->vec[6][2]; 02155 02156 ml_p= MEM_mallocN(sizeof(ml_pointer), "ml_pointer"); 02157 ml_p->ml= mainb[a]; 02158 BLI_addtail(&node->elems, ml_p); 02159 02160 if(mainb[a]->flag & MB_NEGATIVE) { 02161 /* number of negative MetaElem in scene */ 02162 metaball_tree->neg++; 02163 } 02164 else{ 02165 /* number of positive MetaElem in scene */ 02166 metaball_tree->pos++; 02167 } 02168 } 02169 02170 /* size of first node */ 02171 size[0]= node->x_max - node->x_min; 02172 size[1]= node->y_max - node->y_min; 02173 size[2]= node->z_max - node->z_min; 02174 02175 /* first node is subdivided recursively */ 02176 subdivide_metaball_octal_node(node, size[0], size[1], size[2], metaball_tree->depth); 02177 } 02178 02179 void metaball_polygonize(Scene *scene, Object *ob, ListBase *dispbase) 02180 { 02181 PROCESS mbproc; 02182 MetaBall *mb; 02183 DispList *dl; 02184 int a, nr_cubes; 02185 float *ve, *no, totsize, width; 02186 02187 mb= ob->data; 02188 02189 if(totelem==0) return; 02190 if(!(G.rendering) && (mb->flag==MB_UPDATE_NEVER)) return; 02191 if(G.moving && mb->flag==MB_UPDATE_FAST) return; 02192 02193 curindex= totindex= 0; 02194 indices= NULL; 02195 thresh= mb->thresh; 02196 02197 /* total number of MetaElems (totelem) is precomputed in find_basis_mball() function */ 02198 mainb= MEM_mallocN(sizeof(void *)*totelem, "mainb"); 02199 02200 /* initialize all mainb (MetaElems) */ 02201 totsize= init_meta(scene, ob); 02202 02203 if(metaball_tree){ 02204 free_metaball_octal_node(metaball_tree->first); 02205 MEM_freeN(metaball_tree); 02206 metaball_tree= NULL; 02207 } 02208 02209 /* if scene includes more then one MetaElem, then octal tree optimalisation is used */ 02210 if((totelem > 1) && (totelem <= 64)) init_metaball_octal_tree(1); 02211 if((totelem > 64) && (totelem <= 128)) init_metaball_octal_tree(2); 02212 if((totelem > 128) && (totelem <= 512)) init_metaball_octal_tree(3); 02213 if((totelem > 512) && (totelem <= 1024)) init_metaball_octal_tree(4); 02214 if(totelem > 1024) init_metaball_octal_tree(5); 02215 02216 /* don't polygonize metaballs with too high resolution (base mball to small) 02217 * note: Eps was 0.0001f but this was giving problems for blood animation for durian, using 0.00001f */ 02218 if(metaball_tree) { 02219 if( ob->size[0] <= 0.00001f * (metaball_tree->first->x_max - metaball_tree->first->x_min) || 02220 ob->size[1] <= 0.00001f * (metaball_tree->first->y_max - metaball_tree->first->y_min) || 02221 ob->size[2] <= 0.00001f * (metaball_tree->first->z_max - metaball_tree->first->z_min)) 02222 { 02223 new_pgn_element(-1); /* free values created by init_meta */ 02224 02225 MEM_freeN(mainb); 02226 02227 /* free tree */ 02228 free_metaball_octal_node(metaball_tree->first); 02229 MEM_freeN(metaball_tree); 02230 metaball_tree= NULL; 02231 02232 return; 02233 } 02234 } 02235 02236 /* width is size per polygonize cube */ 02237 if(G.rendering) width= mb->rendersize; 02238 else { 02239 width= mb->wiresize; 02240 if(G.moving && mb->flag==MB_UPDATE_HALFRES) width*= 2; 02241 } 02242 /* nr_cubes is just for safety, minimum is totsize */ 02243 nr_cubes= (int)(0.5f+totsize/width); 02244 02245 /* init process */ 02246 mbproc.function = metaball; 02247 mbproc.size = width; 02248 mbproc.bounds = nr_cubes; 02249 mbproc.cubes= NULL; 02250 mbproc.delta = width/(float)(RES*RES); 02251 02252 polygonize(&mbproc, mb); 02253 02254 MEM_freeN(mainb); 02255 02256 /* free octal tree */ 02257 if(totelem > 1){ 02258 free_metaball_octal_node(metaball_tree->first); 02259 MEM_freeN(metaball_tree); 02260 metaball_tree= NULL; 02261 } 02262 02263 if(curindex) { 02264 dl= MEM_callocN(sizeof(DispList), "mbaldisp"); 02265 BLI_addtail(dispbase, dl); 02266 dl->type= DL_INDEX4; 02267 dl->nr= mbproc.vertices.count; 02268 dl->parts= curindex; 02269 02270 dl->index= indices; 02271 indices= NULL; 02272 02273 a= mbproc.vertices.count; 02274 dl->verts= ve= MEM_mallocN(sizeof(float)*3*a, "mballverts"); 02275 dl->nors= no= MEM_mallocN(sizeof(float)*3*a, "mballnors"); 02276 02277 for(a=0; a<mbproc.vertices.count; a++, no+=3, ve+=3) { 02278 ve[0]= mbproc.vertices.ptr[a].position.x; 02279 ve[1]= mbproc.vertices.ptr[a].position.y; 02280 ve[2]= mbproc.vertices.ptr[a].position.z; 02281 02282 no[0]= mbproc.vertices.ptr[a].normal.x; 02283 no[1]= mbproc.vertices.ptr[a].normal.y; 02284 no[2]= mbproc.vertices.ptr[a].normal.z; 02285 } 02286 } 02287 02288 freepolygonize(&mbproc); 02289 } 02290