/* Copyright (C) 2012 University of Oxford */ /* Part of FSL - FMRIB's Software Library http://www.fmrib.ox.ac.uk/fsl fsl@fmrib.ox.ac.uk Developed at FMRIB (Oxford Centre for Functional Magnetic Resonance Imaging of the Brain), Department of Clinical Neurology, Oxford University, Oxford, UK LICENCE FMRIB Software Library, Release 5.0 (c) 2012, The University of Oxford (the "Software") The Software remains the property of the University of Oxford ("the University"). The Software is distributed "AS IS" under this Licence solely for non-commercial use in the hope that it will be useful, but in order that the University as a charitable foundation protects its assets for the benefit of its educational and research purposes, the University makes clear that no condition is made or to be implied, nor is any warranty given or to be implied, as to the accuracy of the Software, or that it will be suitable for any particular purpose or for use under any specific conditions. Furthermore, the University disclaims all responsibility for the use which is made of the Software. It further disclaims any liability for the outcomes arising from using the Software. The Licensee agrees to indemnify the University and hold the University harmless from and against any and all claims, damages and liabilities asserted by third parties (including claims for negligence) which arise directly or indirectly from the use of the Software or the sale of any products based on the Software. No part of the Software may be reproduced, modified, transmitted or transferred in any form or by any means, electronic or mechanical, without the express permission of the University. The permission of the University is not required if the said reproduction, modification, transmission or transference is done without financial return, the conditions of this Licence are imposed upon the receiver of the product, and all original and amended source code is included in any transmitted product. You may be held legally responsible for any copyright infringement that is caused or encouraged by your failure to abide by these terms and conditions. You are not permitted under this Licence to use this Software commercially. Use for which any financial return is received shall be defined as commercial use, and includes (1) integration of all or part of the source code or the Software into a product for sale or license by or on behalf of Licensee to third parties or (2) use of the Software or any derivative of it for research with the final aim of developing software products for sale or license to a third party or (3) use of the Software or any derivative of it for research with the final aim of developing non-software products for sale or license to a third party, or (4) use of the Software to provide any service to an external organisation for which payment is received. If you are interested in using the Software commercially, please contact Isis Innovation Limited ("Isis"), the technology transfer company of the University, to negotiate a licence. Contact details are: innovation@isis.ox.ac.uk quoting reference DE/9564. */ #include #include #include //STL includes #include #include #include using namespace std; #define PI 3.141592653589793238462643383 //using namespace fslsurface_name; namespace fslsurface_name { /* void normalize( float3 & v1) { float mag=sqrtf(v1.x*v1.x + v1.y*v1.y + v1.z*v1.z); v1.x/=mag; v1.y/=mag; v1.z/=mag; } */ float dot(const vec3 & v1, const vec3 & v2 ) { return (v1.x*v2.x + v1.y*v2.y + v1.z*v2.z); } vec3 cross(const vec3 & v1, const vec3 & v2 ) { return vec3( v1.x*v2.y - v1.y*v2.x , \ v1.y*v2.z - v1.z*v2.y , \ v1.z*v2.x - v1.x*v2.z); } vec3 subtract(const vec3 & v1, const vec3 & v2) { return vec3(v1.x - v2.x, v1.y - v2.y, v1.z - v2.z ); } template void projectVectorsOntoNormals(fslSurface& surf ,const unsigned int & index) { //caculate normals and normalize surf.calculateNormals(true,false); vector new_vecs(surf.N_vertices*3); vector new_scalars(surf.N_vertices); string name= surf.vector_names[index]; surf.vector_names.push_back(name + "_proj_onto_normals"); typename vector::iterator i_vec = surf.vector_data[index].begin(); typename vector::iterator i_vec_new = new_vecs.begin(); typename vector::iterator i_sc_new = new_scalars.begin(); for ( typename vector< vertex >::iterator i_v = surf.vbegin(); i_v != surf.vend(); ++i_v, i_vec+=3, i_vec_new+=3,++i_sc_new) { T mag = i_v->nx * (*i_vec) + i_v->ny * (*(i_vec+1)) + i_v->nz * (*(i_vec+2)); *i_sc_new = mag; *i_vec_new = i_v->nx * mag; *(i_vec_new+1) = i_v->ny * mag; *(i_vec_new+2) = i_v->nz * mag; } surf.scalar_data.insert(surf.scalar_data.begin(),new_scalars); surf.scalar_names.insert(surf.scalar_names.begin(), name + "_proj_onto_normals"); surf.vector_data.insert(surf.vector_data.begin(),new_vecs); surf.vector_names.insert(surf.vector_names.begin(), name + "_proj_onto_normals"); } template void projectVectorsOntoNormals(fslSurface& surf ,const unsigned int & index); template void apply_xfm( fslSurface& surf, const string & xfm ){ vector vxfm(16,0); ifstream fmat(xfm.c_str()); T temp; for (typename vector::iterator i_v = vxfm.begin(); i_v != vxfm.end(); ++i_v) { fmat>>temp; *i_v = temp; } fmat.close(); apply_xfm(surf,vxfm ); } template void apply_xfm( fslSurface& surf, const string & xfm ); template void apply_xfm( fslSurface& surf, const string & xfm ); template void apply_xfm( fslSurface& surf, const string & xfm ); template void apply_xfm( fslSurface& surf, const string & xfm ); template void apply_xfm( fslSurface& surf, const vector & xfm ) { //assume bottom row of xfm is 0 0 0 1 for ( typename vector< vertex >::iterator i_vert = surf.vertices.begin(); i_vert != surf.vertices.end(); ++i_vert ) { T x_in = i_vert->x; T y_in = i_vert->y; T z_in = i_vert->z; i_vert->x = xfm[0] * x_in + xfm[1] * y_in + xfm[2] * z_in + xfm[3] ; i_vert->y = xfm[4] * x_in + xfm[5] * y_in + xfm[6] * z_in + xfm[7] ; i_vert->z = xfm[8] * x_in + xfm[9] * y_in + xfm[10] * z_in + xfm[11] ; } } template void apply_xfm( fslSurface& surf, const vector & xfm ); template void apply_xfm( fslSurface& surf, const vector & xfm ); template void apply_xfm( fslSurface& surf, const vector & xfm ); template void apply_xfm( fslSurface& surf, const vector & xfm ); template void apply_xfm( fslSurface& surf, const vector & xfm ); template void apply_warp( fslSurface& surf, const NEWIMAGE::volume4D & warp ) { vec3 dims= vec3(warp.xdim(),warp.ydim(),warp.zdim()); for ( typename vector< vertex >::iterator i_v=surf.vbegin(); i_v != surf.vend(); ++i_v) { // cout<x<<" "<y<<" "<z<<" "<< vec3 vox_coord = vec3(i_v->x / dims.x,i_v->y / dims.y, i_v->z / dims.z ); i_v->x += warp[0].interpolate(vox_coord.x,vox_coord.y,vox_coord.z); i_v->y += warp[1].interpolate(vox_coord.x,vox_coord.y,vox_coord.z); i_v->z += warp[2].interpolate(vox_coord.x,vox_coord.y,vox_coord.z); } } template void apply_warp( fslSurface& surf, const NEWIMAGE::volume4D & warp ); template void deformSurface( fslSurface& surf , const unsigned int & N, const T & w_norm, const T & w_tan ) { surf.calculateNormals(); if (surf.adj_tris.empty()) surf.computeAdjacentTriangles(); if (surf.adj_verts.empty()) surf.computeAdjacentVertices(true); for ( unsigned int iter = 0 ; iter < N ; ++iter ) { for ( typename vector< vertex >::iterator i_vert = surf.vertices.begin() ; i_vert != surf.vertices.end(); ++i_vert ) { i_vert->x += i_vert->nx * w_norm; i_vert->y += i_vert->ny * w_norm; i_vert->z += i_vert->nz * w_norm; } } } template void remove_connection( vector< list > & connections, const unsigned int & index ) { //remove backward connections for (typename list::iterator i = connections[index].begin() ; i != connections[index].end(); ++i) for (typename list::iterator ii = connections[*i].begin() ; ii != connections[*i].end(); ++ii) if (*i == index) { connections[index].erase(ii); break; } //clear row (forward conections) connections[index].clear(); } template void remove_item( list & connections, const unsigned int & index ) { //remove backward connections for (typename list::iterator i = connections.begin() ; i != connections.end(); ++i) if (*i == index) { connections.erase(i); break; } } //this computes the furhter down connections template vector< list > computeNextLevelConnections( vector< std::list > & connections0, vector< std::list > & connections1 ) { vector< list > new_cons(connections0.size()); unsigned int curr_ind = 0 ; typename vector< list >::iterator i_new = new_cons.begin(); for ( typename vector< list >::iterator i0 = connections0.begin(); i0 != connections0.end(); ++i0,++i_new,++curr_ind) for ( typename list::iterator ii0 = i0->begin(); ii0 != i0->end(); ++ii0) { for ( typename list::iterator ii1 = connections1[*ii0].begin(); ii1 != connections1[*ii0].end(); ++ii1) { if (curr_ind != *ii1) { //check too make sure it does not map back onto a connection in the first matrix bool exists=false; for ( typename list::iterator iii0 = i0->begin(); iii0 != i0->end(); ++iii0) if(*iii0 == *ii1) exists=true; //make sure it does already exist in two double connection for ( typename list::iterator ii_new = i_new->begin(); ii_new != i_new->end(); ++ii_new) if(*ii_new == *ii1) exists=true; if (!exists) i_new->push_back(*ii1); } } } // cout<<"new cons"< >::iterator i0 = new_cons.begin(); i0 != new_cons.end(); ++i0) // { // for ( typename list::iterator ii0 = i0->begin(); ii0 != i0->end(); ++ii0) // cout<<*ii0<<" "; // cout< > computeNextLevelConnections( vector< std::list > & connections0, vector< std::list > & connections1 ); //this computes the furhter down connections template vector< list< pair > > computeNextLevelConnectionsWithMultiWayPoints( vector< std::list > & connections0, vector< std::list > & connections1 ) { //second in pair vector< list< pair > > new_cons(connections0.size()); // vector< list< pair > > waypoints(); unsigned int curr_ind = 0 ; typename vector< list< pair > >::iterator i_new = new_cons.begin(); for ( typename vector< list >::iterator i0 = connections0.begin(); i0 != connections0.end(); ++i0,++i_new,++curr_ind) for ( typename list::iterator ii0 = i0->begin(); ii0 != i0->end(); ++ii0) { for ( typename list::iterator ii1 = connections1[*ii0].begin(); ii1 != connections1[*ii0].end(); ++ii1) { if (curr_ind != *ii1) { //check too make sure it does not map back onto a connection in the first matrix bool exists=false; for ( typename list::iterator iii0 = i0->begin(); iii0 != i0->end(); ++iii0) if(*iii0 == *ii1) exists=true; //make sure it does already exist in two double connection // for ( typename list< pair >::iterator ii_new = i_new->begin(); ii_new != i_new->end(); ++ii_new) // if(ii_new->first == *ii1) // exists=true; if (!exists) i_new->push_back( pair(*ii1,*ii0)); } } } // cout<<"new cons"< > >::iterator i0 = new_cons.begin(); i0 != new_cons.end(); ++i0, ++c) // { // for ( typename list< pair >::iterator ii0 = i0->begin(); ii0 != i0->end(); ++ii0) // cout<<"conn2 "<first<<" waypoint: "<second< > > computeNextLevelConnectionsWithMultiWayPoints( vector< std::list > & connections0, vector< std::list > & connections1 ); template void cluster( fslSurface& surf , const unsigned int & index, const T & threshold ) { if (surf.adj_verts_bi) surf.adj_verts.clear(); if (surf.adj_verts.empty()) surf.computeAdjacentVertices(false); //copy over adjacent vertices //going to delete as search through scalars vector clusters(surf.N_vertices,0); vector used(surf.N_vertices,0); //keep track of which cluster were on T cluster_count=1; typename vector::iterator i_cl = clusters.begin(); vector< list > connections = surf.adj_verts; unsigned int Nconns = connections.size(); if (surf.N_vertices != Nconns) throw fslSurfaceException( "cluster : Number of connection does not match number of vertices" ); typename vector::const_iterator i_sc = surf.const_scbegin(index); // for ( typename vector::const_iterator i_sc = surf.const_scbegin(index); i_sc != surf.const_scbegin(index); ++i_sc, ++i_cl, ++i_con, ++v_index) for ( T2 v_index = 0 ; v_index < Nconns; ++v_index) { // cout<<"v_index "<= threshold) //will enter this once per cluster { // clusters[v_index]=cluster_count; list stack; stack.push_back(v_index); while (!stack.empty()){ T2 current_vert = stack.back(); //cout<<"stack "<::const_iterator i_con = connections[current_vert].begin(); i_con != connections[current_vert].end(); ++i_con) { //if scalar is above threshold //cout<<"val "<<*i_con<<" "<<*(i_sc + *i_con)<<" "<= threshold) { stack.push_back(*i_con); remove_item(connections[*i_con], current_vert); } else//other remove connection with the cluster { fslsurface_name::remove_connection(connections, *i_con); // used[*i_con]=1; } } clusters[current_vert] = cluster_count; //vert has beeen set and connections pushed to stack remove_connection(connections, current_vert); // used[current_vert]; } ++cluster_count; }else { //delete all connections fslsurface_name::remove_connection(connections,v_index); used[v_index]; } } } // surf.scalar_data.insert(surf.scalar_data.begin(),clusters); surf.scalar_data.push_back(clusters); stringstream ss; ss<>thresh; surf.scalar_names.push_back("clusters_"+thresh); } template void cluster( fslSurface& surf , const unsigned int & index, const float & threshold ); template void sc_smooth_mean_neighbour( fslSurface& surf , const unsigned int & index ) { cout<<"average using mean neighbours"< sc_mean = surf.scalar_data[index] ; //(surf.getNumberOfVertices(),0); vector< list > connections = surf.adj_verts; typename vector< list >::iterator i_adj = connections.begin(); typename std::vector< T >::const_iterator i_sc_orig = surf.const_scbegin(index); // vector::iterator i_scN = sc_mean_N.begin(); for ( vector::iterator i_sc = sc_mean.begin(); i_sc != sc_mean.end(); ++i_sc , ++i_adj) { for (typename list::iterator ii_adj = i_adj->begin(); ii_adj != i_adj->end(); ++ii_adj) { *i_sc += *(i_sc_orig+*ii_adj); // cout<<*(i_sc_orig+*ii_adj)<size()+1; } surf.scalar_data.insert( surf.scalar_data.begin(),sc_mean); surf.scalar_names.insert( surf.scalar_names.begin(),surf.scalar_names[index] + "_mean_neighbour"); } template void sc_smooth_mean_neighbour( fslSurface& surf , const unsigned int & index ); //return the connection plus the distance template //std::list< pair > void conn2dist(fslSurface& surf,const T2 & vert_index, std::list< std::pair > & conns, vector< std::list< pair > > & index_distances ) { //first is connection //second is waypoint //assumes that we traverse the surface from 0 -> N-1 typename vector< vertex >::const_iterator v_iter = surf.const_vbegin(); // cout<<"COMPUTE DIST "< > distances; //look for duplicate conns while (! conns.empty()) { //cout<<"loop"< point = conns.front(); // cout< >::iterator i = index_distances//[vert_index].begin(); i != index_distances[vert_index].end(); ++i) // if ( i->first == point.first ) // {//check to see if exists // not_already_computed = false; // break; // } // if (not_already_computed) // cout<<"pointnew "<N-i1 if (point.first > vert_index) { bool notFoundDup=true; float distance = surf.L2norm(vert_index,point.second ) + surf.L2norm(point.second,point.first ); for (typename list< std::pair >::iterator i = conns.begin(); i != conns.end();++i) { // cout<<"found duplicate???? "<first<<" "<second<first ) { // cout<<"found duplicate "<first<<" "<second< rot_axis = surf.subtractVerts(i->second, point.second); // normalize(rot_axis); rot_axis.normalize(); //subtraction of both way points is // cout<<"rot axis "< P0 = surf.getVertexCoord(vert_index); vec3 P1 = surf.getVertexCoord(point.second); vec3 P2 = surf.getVertexCoord(point.first); vec3 P3 = surf.getVertexCoord(i->second); float p10_dot_rot = ((P0.x-P1.x)*rot_axis.x + (P0.y-P1.y)*rot_axis.y + (P0.z-P1.z)*rot_axis.z); //projection // cout<<"P01do Rot "< base_perp = vec3(P1.x + rot_axis.x * p10_dot_rot - P0.x , \ P1.y + rot_axis.y * p10_dot_rot - P0.y , \ P1.z + rot_axis.z * p10_dot_rot - P0.z ); // cout<<"base_perp "< p12 = subtract(P2, P1); float p12_dot_rot = dot(p12,rot_axis); // float p12_dot_rot = ((P2.x-P1.x)*rot_axis.x + (P2.y-P1.y)*rot_axis.y + (P2.z-P1.z)*rot_axis.z); float p12_dot_perp = dot(p12,base_perp); vec3 P2_rotated = vec3(P1.x + rot_axis.x * p12_dot_rot + base_perp.x * p12_dot_perp, \ P1.y + rot_axis.y * p12_dot_rot + base_perp.y * p12_dot_perp, \ P1.z + rot_axis.z * p12_dot_rot + base_perp.z * p12_dot_perp); ///------------------TEST FOR INTERSECTION WITH ROTATION EDGES---------------// //do a copplanar check all points should be in same plane, don't use rotated point in case //numerical error. Use a way point and the original, one from each line float t1 =0.0f; //make sure its not in the Z plane (i.e. x or y are all the same) //simplify, dont need to test P2, because will by definition be inplabe with 0,1,2 //usinf roateted could have introduced numerical errors if ((P0.x==P1.x)&&(P1.x==P3.x)) t1 = ((P1.z - P0.z)*(P2_rotated.y-P0.y) - (P2_rotated.z-P0.z)*(P1.y-P0.y))/( (P2_rotated.z-P0.z)*(P3.y-P1.y) - (P3.z-P1.z)*(P2_rotated.y-P0.y) ); else if ((P0.y==P1.y)&&(P1.y==P3.y)) t1 = ((P1.x - P0.x)*(P2_rotated.z-P0.z) - (P2_rotated.x-P0.x)*(P1.z-P0.z))/( (P2_rotated.x-P0.x)*(P3.z-P1.z) - (P3.x-P1.x)*(P2_rotated.z-P0.z) ); else t1 = ((P1.x - P0.x)*(P2_rotated.y-P0.y) - (P2_rotated.x-P0.x)*(P1.y-P0.y))/( (P2_rotated.x-P0.x)*(P3.y-P1.y) - (P3.x-P1.x)*(P2_rotated.y-P0.y) ); //*************DEBUGGING*******************// // cout<<"set roatated point "<second<<" "<=0)&& (t1<=1)) { // cout<<"Protated "<(vert_index,distance)); // distances. index_distances[vert_index].push_back( pair(point.first, distance ) ); // } } } // return distances; } template void sc_smooth_gaussian_geodesic( fslSurface& surf ,const unsigned int & index, const T & sigma , const T & extent , bool run4D ) { cout<<"do smooth"< > distance_1; surf.adj_verts.clear(); if (surf.adj_verts.empty()) surf.computeAdjacentVertices(true);//only use double sided, calculate distance in each direction //for testing start with onnection is dist 1 //start with just single neighbours // vector< list > connections = surf.adj_verts; // vector< list > connections_2 = computeNextLevelConnectionsWithMultiWayPoints(surf.adj_verts,surf.adj_verts); vector< list< pair > > connections_2 = computeNextLevelConnectionsWithMultiWayPoints(surf.adj_verts,surf.adj_verts); vector< list< pair > > index_distance(surf.N_vertices); // bool not_already_computed=true; // for ( typename std::list< pair >::iterator i = index_distances[vert_index].begin(); i != index_distances[vert_index].end(); ++i) // if ( i->first == point.first ) // {//check to see if exists // not_already_computed = false; // break; // } //if (not_already_computed) //{ {//do 1-neighbour ring unsigned int curr_index=0; typename vector< list< pair > >::iterator i_id = index_distance.begin(); for ( typename vector< list >::iterator i_adj = surf.adj_verts.begin(); i_adj != surf.adj_verts.end();++i_adj,++curr_index,++i_id) { // cout<<"current index "<::iterator ii_adj = i_adj->begin(); ii_adj != i_adj->end(); ++ii_adj) { if (curr_index<*ii_adj) { float dist = surf.L2norm(curr_index,*ii_adj); // cout<<"hmm "<size()<<" "<push_back( pair(*ii_adj, dist ) ); index_distance[*ii_adj].push_back( pair(curr_index,dist) ); } } // cout< > >::iterator i_adj = index_distance.begin(); i_adj != index_distance.end();++i_adj ) // { // for (typename list< pair >::iterator ii_adj = i_adj->begin(); ii_adj != i_adj->end(); ++ii_adj) // { // cout<first<<","<second<<" "; // } // cout< > >::iterator i_id = index_distance.begin(); for ( typename vector< list< pair > >::iterator i_adj = connections_2.begin(); i_adj != connections_2.end();++i_adj,++curr_index,++i_id) { // cout<<"current index21 "<size()<push_back( pair(1,1)); //list< pair > con2s = conn2dist(surf, curr_index, *i_adj, index_distance); //i_id->merge(con2s); //break; // cout<<"current index2 "< >::iterator ii_adj = i_adj->begin(); ii_adj != i_adj->end(); ++ii_adj) // { // cout<<"hmm2 "<size()<first<<" "<second<<" "<first) <push_back( pair(ii_adj->first, surf.conn2dist(curr_index,ii_adj->first) ) ); // } // cout< > >::iterator i_adj = index_distance.begin(); i_adj != index_distance.end();++i_adj ) // { // for (typename list< pair >::iterator ii_adj = i_adj->begin(); ii_adj != i_adj->end(); ++ii_adj) // { // cout<first<<","<second<<" "; /// } // cout< graph;//(index_distance.size(),1e16);//the distance vector< bool > graph_visited;//(index_distance.size(),false);//the distance vector< bool > in_queue;//(index_distance.size(),false);//the distance //do firts vertex then will follow with the rest // unsigned int curr_index=0; //for a given node visit all negihbours //this will be in large loop around all vertices // T2 vertex_index = 0; T2 Nverts = index_distance.size(); // Nverts=1; vector< vector > smoothed_scalars(1,vector(Nverts,0)); vector smoothed_names(surf.getNumberOfScalarData()); // vector< vector > smoothed_scalars4D; if (run4D) { smoothed_scalars.assign(surf.getNumberOfScalarData(),vector(Nverts,0)); }else{ } //first smoothed vectors typename vector::iterator i_sc_smoothed = smoothed_scalars[0].begin(); for ( T2 vertex_index =0 ; vertex_index < Nverts; ++vertex_index, ++i_sc_smoothed) { if ( vertex_index%100 == 0 ) cout<<"Smoothing at vertex "<::iterator i = graph.begin(); i != graph.end(); ++i,++count2) // cout<::iterator i_graph = graph.begin()+vertex_index; list< pair > points_to_visit;//with distance typename list< pair >::iterator i_neigh_end = index_distance[vertex_index].end(); //this is actually just copying info from center (directtly from index_distance) //this is for immediate neighbours for (typename list< pair >::iterator i_neigh = index_distance[vertex_index].begin(); i_neigh != index_distance[vertex_index].end(); ++i_neigh) { graph[i_neigh->first] = i_neigh->second;//dist; points_to_visit.push_back( pair(i_neigh->second,i_neigh->first) ); } graph_visited[vertex_index] = true; // cout<<"graph graph graph "<::iterator i = graph.begin(); i != graph.end(); ++i,++count2) // cout< p_curr = points_to_visit.front(); // cout<<"new point of queue "< >::iterator i_neigh_end = index_distance[curr_index].end(); //this is actually just copying info from center (directtly from index_distance) for (typename list< pair >::iterator i_neigh = index_distance[curr_index].begin(); i_neigh != i_neigh_end; ++i_neigh) { //set current distance of node to new distance plus distance at current node //don't need to visit if its all ready been totally added if (!graph_visited[i_neigh->first]) { // cout<<"enter if "<first<second; if (graph[i_neigh->first] > dist) { // cout<<"set distance "<second<<" "<first<<" "<first]<<" "<first] = dist; } if (dist<= extent) { //if the neighbour has not been marked as visited if ((!in_queue[i_neigh->first])) { points_to_visit.push_back( pair (dist,i_neigh->first) ); in_queue[i_neigh->first]=true; }else { //update distance of point in queue for (typename list< pair >::iterator i = points_to_visit.begin(); i!= points_to_visit.end(); ++i) { if ( i->second == curr_index) { i->first = dist; break; } } } } } } graph_visited[curr_index] = true; } //cout<<"graph results "<::iterator i = graph.begin(); i != graph.end(); ++i,++count) // cout< >::iterator ii_sc = surf.scalar_data.begin(); ii_sc != surf.scalar_data.end(); ++ii_sc,++count) { if (vertex_index==0) { stringstream ss; string ssigma,sextent; ss<>ssigma; ss<>sextent; // cout<<"setname "<<(surf.getScalarName(count)+"_smoothed_sigma"+ssigma+"_extent"+sextent)<::iterator i_sc = ii_sc->begin(); for ( typename vector::iterator i_g = graph.begin(); i_g != graph.end(); ++i_g,++i_sc ) { //make sure were within the extent if ( *i_g <= extent ) { float factor = expf( -0.5*(*i_g)*(*i_g)*invsigma2 ); *i_sc_smoothed += (*i_sc)*factor; norm+=factor; } } *i_sc_smoothed/=norm; } }else{ float norm=0; typename vector::iterator i_sc = surf.scalar_data[index].begin(); for ( typename vector::iterator i_g = graph.begin(); i_g != graph.end(); ++i_g,++i_sc ) { //make sure were within the extent if ( *i_g <= extent ) { float factor = expf( -0.5*(*i_g)*(*i_g)*invsigma2 ); *i_sc_smoothed += (*i_sc)*factor; norm+=factor; } } *i_sc_smoothed/=norm; } // cout<<"norm "<>ssigma; ss<>sextent; if (run4D) { // cout<<"replace scalars "< >::iterator i_adj = connections.begin(); //while ( !i_adj->empty() ) { // } // for (typename list::iterator ii_adj = i_adj->begin(); ii_adj != i_adj->end(); ++ii_adj) // { // *i_sc += *(i_sc_orig+*ii_adj); // // cout<<*(i_sc_orig+*ii_adj)<& surf ,const unsigned int & index, const float & variance , const float & extent , bool run4D); template void multiplyScalarsByScalars(fslSurface& surf , const unsigned int & index, const fslSurface& surf2 , const unsigned int & index1) { vector new_sc( surf.N_vertices , 0 ); typename vector::iterator i_sc_new= new_sc.begin(); typename vector::const_iterator i_sc1= surf2.const_scbegin(index); for ( typename vector::const_iterator i_sc0= surf.const_scbegin(index); i_sc0 != surf.const_scend(index); ++i_sc0, ++i_sc1, ++i_sc_new) { *i_sc_new = (*i_sc0) * (*i_sc1); } surf.insertScalars(new_sc, 0, "sc_mul_sc"); } template void multiplyScalarsByScalars( fslSurface& surf ,const unsigned int & index, const fslSurface& surf2 ,const unsigned int & index1 ); template void subtractScalarsFromScalars(fslSurface& surf , const unsigned int & index, const fslSurface& surf2 , const unsigned int & index1, string name ){ vector new_sc( surf.N_vertices , 0 ); typename vector::iterator i_sc_new= new_sc.begin(); typename vector::const_iterator i_sc1= surf2.const_scbegin(index); for ( typename vector::const_iterator i_sc0= surf.const_scbegin(index); i_sc0 != surf.const_scend(index); ++i_sc0, ++i_sc1, ++i_sc_new) { cout<<"sc "<<*i_sc0<<" "<<*i_sc1<& surf ,const unsigned int & index, const fslSurface & surf2 ,const unsigned int & index1,string name ); template void runMarchingCubesOnAllLabels( fslSurface& surf, const T3* image, const fslsurface_name::image_dims & dims, const T3 & thresh_in) { set all_labels; int nvoxels = static_cast(dims.xsize * dims.ysize * dims.zsize); for (int i = 0 ; i < nvoxels; ++i) all_labels.insert(image[i]); for ( typename set::iterator i_lb = all_labels.begin(); i_lb != all_labels.end(); ++i_lb) if (*i_lb) marchingCubes(surf, image, dims, static_cast (*i_lb),*i_lb,EQUAL_TO); } template void runMarchingCubesOnAllLabels( fslSurface& surf, const float* image, const fslsurface_name::image_dims & dims, const float & thresh_in); template void runMarchingCubesOnAllLabels( fslSurface& surf, const int* image, const fslsurface_name::image_dims & dims, const int & thresh_in); template void marchingCubes( fslSurface& surf, const T3* image, const fslsurface_name::image_dims & dims, const T3 & thresh_in, const T & label, const MarchingCubesMode & mode) { // cout<<"Marching cubes input "< im(3,3,3); int tc=0; for (int z=0;z<3;++z) for (int y=0;y<3;++y) for (int x=0;x<3;++x,++tc) im.value(x,y,z)=tc; dims.xsize=3; dims.ysize=3; dims.zsize=3; */ //-------------TEST CODE--------------// // vertices.clear(); // faces.clear(); // scalar_data.clear(); // vector_data.clear(); if (surf.scalar_data.empty()) { surf.scalar_names.push_back("mcubes_labels"); vector sc_temp; surf.scalar_data.push_back(sc_temp); } int ystride = static_cast(dims.xsize); int zstride = static_cast(dims.xsize * dims.ysize); int nvoxels = static_cast(dims.xsize * dims.ysize * dims.zsize); T3 thresh = thresh_in; //cout<<"apply threshold "<(image[i]) == static_cast(thresh) ) ? 1 : 0)<(image[i]) == static_cast(thresh) ) ? 1 : 0) == thresh ) // cout<<"image match "<(image[i]+0.5) == static_cast(thresh+0.5) ) ? 1 : 0 ; } thresh=0.5; // imth[i] = ( image[i] == thresh ) ? image[i] : 0 ; }else if (mode == GREATER_THAN) { for (int i = 0 ; i < nvoxels; ++i) imth[i] = ( image[i] > thresh ) ? image[i] : 0 ; }else if (mode == GREATER_THAN_EQUAL_TO) { for (int i = 0 ; i < nvoxels; ++i) imth[i] = ( image[i] >= thresh ) ? image[i] : 0 ; } //cout<<"done apply threshold "< thresh ) ? im_ptr[i] : 0 ; */ //-------------TEST CODE--------------// // //cout<<"Done thresholding the image "< > index_2_index_all; for ( int z = 0 ; z < dims.zsize-1; ++z) { //cout<<"xyz "<(x + y*ystride + z*zstride)<<". "<(x + y*ystride + z*zstride)]<(x + y*ystride + z*zstride)]<(x + y*ystride + z*zstride)] cubeindex = 0; /* vertVals[0] = imth[x + (y+1)*ystride + z*zstride ]; vertVals[1] = imth[(x+1) + (y+1)*ystride + z*zstride ]; vertVals[2] = imth[(x+1) + y*ystride + z*zstride ]; vertVals[3] = imth[x + (y+1)*ystride + z*zstride ]; vertVals[4] = imth[x + (y+1)*ystride + z*zstride ]; vertVals[5] = imth[x + (y+1)*ystride + z*zstride ]; vertVals[6] = imth[x + (y+1)*ystride + z*zstride ]; vertVals[7] = imth[x + (y+1)*ystride + z*zstride ]; */ vertVals[0] = imth[index + ystride ]; vertVals[1] = imth[index + ystride +1]; vertVals[2] = imth[index + 1]; vertVals[3] = imth[index ]; vertVals[4] = imth[index +zstride + ystride ]; vertVals[5] = imth[index +zstride + ystride +1]; vertVals[6] = imth[index +zstride + 1 ]; vertVals[7] = imth[index +zstride ]; // if ((z==158)&&(y==153)) // //cout<<"cube vals:"< index_2_index; // if ((z==158)&&(y==153)) // //cout<<"check edge table "<(x + y*ystride + z*zstride)<<". "<(x + y*ystride + z*zstride)]< vnew(label); vertexInterp(thresh, \ vec3(x*dims.xdim,(y+1)*dims.ydim, z*dims.zdim),\ vertVals[0], \ vec3((x+1)*dims.xdim,(y+1)*dims.ydim, z*dims.zdim), \ vertVals[1], vnew ); surf.vertices.push_back(vnew); surf.scalar_data[0].push_back(label); index_2_index.insert( pair (0,vert_index+vert_count) ); ++vert_count; }else { index_2_index.insert( pair(0, index_2_index_all[x +y*ystride + (z-1)*zstride].find(4)->second) ); } } if (edgeTable[cubeindex] & 2) { //edge p1-p2 if (z==0) { vertex vnew(label); vertexInterp(thresh, \ vec3((x+1)*dims.xdim,(y+1)*dims.ydim, z*dims.zdim), \ vertVals[1], \ vec3((x+1)*dims.xdim,y*dims.ydim, z*dims.zdim), \ vertVals[2], vnew ); surf.vertices.push_back(vnew); surf.scalar_data[0].push_back(label); index_2_index.insert( pair (1,vert_index+vert_count) ); ++vert_count; }else { index_2_index.insert( pair(1, index_2_index_all[x +y*ystride + (z-1)*zstride].find(5)->second) ); } }if (edgeTable[cubeindex] & 4) { //edge p2-p3 if ((z==0)&&(y==0)) { vertex vnew(label); vertexInterp(thresh, \ vec3((x+1)*dims.xdim,y*dims.ydim, z*dims.zdim), \ vertVals[2], \ vec3(x*dims.xdim,y*dims.ydim, z*dims.zdim), \ vertVals[3], vnew ); surf.vertices.push_back(vnew); surf.scalar_data[0].push_back(label); index_2_index.insert( pair (2,vert_index+vert_count) ); ++vert_count; }else if (z==0){ index_2_index.insert( pair(2, index_2_index_all[x + (y-1)*ystride ].find(0)->second) ); }else{ index_2_index.insert( pair(2, index_2_index_all[x +y*ystride + (z-1)*zstride].find(6)->second) ); } }if (edgeTable[cubeindex] & 8) { //edge p3-p0 if ((x==0)&&(z==0)) { vertex vnew(label); vertexInterp(thresh, \ vec3(x*dims.xdim,y*dims.ydim, z*dims.zdim), \ vertVals[3], \ vec3(x*dims.xdim,(y+1)*dims.ydim, z*dims.zdim), \ vertVals[0], vnew ); surf.vertices.push_back(vnew); surf.scalar_data[0].push_back(label); index_2_index.insert( pair (3,vert_index+vert_count) ); ++vert_count; }else if (z==0) { index_2_index.insert( pair(3, index_2_index_all[x-1 +y*ystride].find(1)->second) ); }else{ index_2_index.insert( pair(3, index_2_index_all[x +y*ystride + (z-1)*zstride].find(7)->second) ); } } if (edgeTable[cubeindex] & 16) { //edge p4-p5 vertex vnew(label); vertexInterp(thresh, \ vec3(x*dims.xdim,(y+1)*dims.ydim, (z+1)*dims.zdim), \ vertVals[4], \ vec3((x+1)*dims.xdim,(y+1)*dims.ydim, (z+1)*dims.zdim), \ vertVals[5], vnew ); surf.vertices.push_back(vnew); surf.scalar_data[0].push_back(label); index_2_index.insert( pair (4,vert_index+vert_count) ); ++vert_count; } if (edgeTable[cubeindex] & 32) { //edge p1-p2 // //cout<<"edge 5 "< vnew(label); vertexInterp(thresh, \ vec3((x+1)*dims.xdim,(y+1)*dims.ydim, (z+1)*dims.zdim), \ vertVals[5], \ vec3((x+1)*dims.xdim,y*dims.ydim, (z+1)*dims.zdim), \ vertVals[6], vnew ); surf.vertices.push_back(vnew); surf.scalar_data[0].push_back(label); index_2_index.insert( pair (5,vert_index+vert_count) ); ++vert_count; }if (edgeTable[cubeindex] & 64) { //edge p2-p3 if (y==0) { vertex vnew(label); vertexInterp(thresh, \ vec3((x+1)*dims.xdim,y*dims.ydim, (z+1)*dims.zdim), \ vertVals[6], \ vec3(x*dims.xdim,y*dims.ydim, (z+1)*dims.zdim), \ vertVals[7], vnew ); surf.vertices.push_back(vnew); surf.scalar_data[0].push_back(label); index_2_index.insert( pair (6,vert_index+vert_count) ); ++vert_count; }else { index_2_index.insert( pair(6, index_2_index_all[x +(y-1)*ystride + z*zstride].find(4)->second) ); } }if (edgeTable[cubeindex] & 128) { //edge p3-p0 ////cout<<"edge 7 "< vnew(label); vertexInterp(thresh, \ vec3(x*dims.xdim,y*dims.ydim, (z+1)*dims.zdim), \ vertVals[7], \ vec3(x*dims.xdim,(y+1)*dims.ydim, (z+1)*dims.zdim), \ vertVals[4], vnew ); surf.vertices.push_back(vnew); surf.scalar_data[0].push_back(label); index_2_index.insert( pair (7,vert_index+vert_count) ); ++vert_count; } else { ////cout<<"find previous index "<(7, index_2_index_all[(x-1) +y*ystride + z*zstride].find(5)->second) ); } ////cout<<"done edge 7 "< vnew(label); vertexInterp(thresh, \ vec3(x*dims.xdim,(y+1)*dims.ydim, z*dims.zdim), \ vertVals[0], \ vec3(x*dims.xdim,(y+1)*dims.ydim, (z+1)*dims.zdim), \ vertVals[4], vnew ); surf.vertices.push_back(vnew); surf.scalar_data[0].push_back(label); index_2_index.insert( pair (8,vert_index+vert_count) ); ++vert_count; }else { index_2_index.insert( pair(8, index_2_index_all[(x-1) +y*ystride + z*zstride].find(9)->second) ); } } if (edgeTable[cubeindex] & 512){ //edge p1-p5 vertex vnew(label); vertexInterp(thresh, \ vec3((x+1)*dims.xdim,(y+1)*dims.ydim, z*dims.zdim), \ vertVals[1], \ vec3((x+1)*dims.xdim,(y+1)*dims.ydim, (z+1)*dims.zdim), \ vertVals[5], vnew ); surf.vertices.push_back(vnew); surf.scalar_data[0].push_back(label); index_2_index.insert( pair (9,vert_index+vert_count) ); ++vert_count; } if (edgeTable[cubeindex] & 1024){ //edge p2-p6 if (y==0) { vertex vnew(label); vertexInterp(thresh, \ vec3((x+1)*dims.xdim,y*dims.ydim, z*dims.zdim), \ vertVals[2], \ vec3((x+1)*dims.xdim,y*dims.ydim, (z+1)*dims.zdim), \ vertVals[6], vnew ); surf.vertices.push_back(vnew); surf.scalar_data[0].push_back(label); index_2_index.insert( pair (10,vert_index+vert_count) ); ++vert_count; }else { index_2_index.insert( pair(10, index_2_index_all[x +(y-1)*ystride + z*zstride].find(9)->second) ); } } if (edgeTable[cubeindex] & 2048){ //edge p3-p7 if ((y==0)&&(x==0)) { vertex vnew(label); vertexInterp(thresh, \ vec3(x*dims.xdim,y*dims.ydim, z*dims.zdim), \ vertVals[3], \ vec3(x*dims.xdim,y*dims.ydim, (z+1)*dims.zdim), \ vertVals[7], vnew ); surf.vertices.push_back(vnew); surf.scalar_data[0].push_back(label); index_2_index.insert( pair (11,vert_index+vert_count) ); ++vert_count; }else if (x==0) { index_2_index.insert( pair(11, index_2_index_all[x +(y-1)*ystride + z*zstride].find(8)->second) ); }else { index_2_index.insert( pair(11, index_2_index_all[(x-1) +y*ystride + z*zstride].find(10)->second) ); } } index_2_index_all.push_back(index_2_index); // //cout<::iterator i=index_2_index.begin(); i != index_2_index.end(); ++i) // //cout<first<<" "<second<second<<" "<second<<" "<second<second<second<second<second); surf.faces.push_back( index_2_index.find(triTable[cubeindex][i_tri+1])->second); surf.faces.push_back( index_2_index.find(triTable[cubeindex][i_tri+2])->second); } // //cout<<"done triangles"< index_2_index; index_2_index_all.push_back(index_2_index); ++index;//because only going to -1 } map index_2_index; for (int i=0;i::iterator i = faces.begin(); i!=faces.end(); ++i) { // cout<<*i<<" "; // } //cout<<"Number of Triangles "<( fslSurface& surf, const float* image, const fslsurface_name::image_dims & dims, const float & thresh, const float & label, const MarchingCubesMode & mode); template void vertexInterp( const float & thresh, const vec3 & v1 , const float & val1, const vec3 & v2, const float & val2, vertex & vout ) { //vertex vout; //cout<<"interp in ("< ("< void mapScalars(fslSurface& surf, const unsigned int & sc_index, const std::map & sc_map , bool setToCurrentSc ) { vector new_scalars( surf.scalar_data[sc_index].size() ); typename vector::iterator i_new = new_scalars.begin(); typename vector::iterator i_end = surf.scalar_data[sc_index].end(); for ( typename vector::iterator i = surf.scalar_data[sc_index].begin(); i != i_end; ++i, ++i_new) { *i_new = sc_map.find(*i)->second; } surf.scalar_data.push_back(new_scalars); if (setToCurrentSc) { surf.setScalars(surf.scalar_data.size()-1); } } template void mapScalars(fslSurface& surf, const unsigned int & sc_index, const std::map & sc_map , bool setToCurrentSc ); template void mapScalars(fslSurface& surf, const unsigned int & sc_index, const std::map & sc_map , bool setToCurrentSc ); template void maskScalars(fslSurface