本文整理汇总了C++中volume类的典型用法代码示例。如果您正苦于以下问题:C++ volume类的具体用法?C++ volume怎么用?C++ volume使用的例子?那么恭喜您, 这里精选的类代码示例或许可以为您提供帮助。
在下文中一共展示了volume类的20个代码示例,这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞,您的评价将有助于我们的系统推荐出更棒的C++代码示例。
示例1: saveSVMFile
// transforms a 4D Volume in a SVM samples file, based on a mask
void saveSVMFile(volume4D <float> &volSamples, volume<float> &mask, const char *outputFileName, float minValue, vector <int > &indexes, vector <int> &classes)
{
FILE *f;
f = fopen(outputFileName, "wt+");
if (f != NULL)
{
int i, t;
for (int h = 0; h < indexes.size(); h++)
{
// picking the right indexes
t = indexes[h] - 1;
i = 0;
// saving the class value. Remember indexes size is different from classes. classes has the same sime of the number of volumes
fprintf(f, "%d ", classes[t]);
for (int z = 0; z < mask.zsize(); z++)
for (int y = 0; y < mask.ysize(); y++)
for (int x = 0; x < mask.xsize(); x++)
if (mask.value(x, y, z) > minValue)
{
i++;
// writing each voxel value in svm format
fprintf(f, "%d:%f ", i, volSamples.value(x, y, z, t));
}
fprintf(f, "\n");
}
fclose(f);
}
}
开发者ID:InstitutoDOr,项目名称:FriendENGINE,代码行数:29,代码来源:svmfuncs.cpp
示例2: tmp
FnirtFileWriter::FnirtFileWriter(const string& fname,
const volume<float>& fieldx,
const volume<float>& fieldy,
const volume<float>& fieldz)
{
volume<float> tmp(fieldx.xsize(),fieldx.ysize(),fieldx.zsize());
tmp.copyproperties(fieldx);
Matrix aff = IdentityMatrix(4);
common_field_construction(fname,tmp,fieldx,fieldy,fieldz,aff);
}
开发者ID:,项目名称:,代码行数:11,代码来源:
示例3: inside_mesh
volume<float> inside_mesh(const volume<float> & image, const Mesh& m)
{
volume<float> res = image;
int xsize = image.xsize();
int ysize = image.ysize();
int zsize = image.zsize();
volume<short> inside = make_mask_from_mesh(image, m);
for (int k=0; k<zsize; k++)
for (int j=0; j<ysize; j++)
for (int i=0; i<xsize; i++)
res.value(i, j, k) = (1-inside.value(i, j, k)) * image.value(i, j, k);
return res;
}
开发者ID:dungttbk,项目名称:sips,代码行数:13,代码来源:bet2.cpp
示例4: make_mask_from_mesh
//calculates the mask from the mesh by spreading an initial point outside the mesh, and stopping it when the mesh is reached.
volume<short> make_mask_from_mesh(const volume<float> & image, const Mesh& m)
{
// cout<<"make_mask_from_mesh begins"<<endl;
double xdim = (double) image.xdim();
double ydim = (double) image.ydim();
double zdim = (double) image.zdim();
volume<short> mask;
copyconvert(image,mask);
int xsize = mask.xsize();
int ysize = mask.ysize();
int zsize = mask.zsize();
mask = 1;
mask = draw_mesh(mask, m);
vector<Pt> current;
current.clear();
Pt c(0., 0., 0.);
for (vector<Mpoint *>::const_iterator it=m._points.begin(); it!=m._points.end(); it++)
c+=(*it)->get_coord();
c*=(1./m._points.size());
c.X/=xdim; c.Y/=ydim; c.Z/=zdim;
current.push_back(c);
while (!current.empty())
{
Pt pc = current.back();
int x, y, z;
x=(int) pc.X; y=(int) pc.Y; z=(int) pc.Z;
//current.remove(pc);
current.pop_back();
mask.value(x, y, z) = 0;
if (0<=x-1 && mask.value(x-1, y, z)==1) current.push_back(Pt(x-1, y, z));
if (0<=y-1 && mask.value(x, y-1, z)==1) current.push_back(Pt(x, y-1, z));
if (0<=z-1 && mask.value(x, y, z-1)==1) current.push_back(Pt(x, y, z-1));
if (xsize>x+1 && mask.value(x+1, y, z)==1) current.push_back(Pt(x+1, y, z));
if (ysize>y+1 && mask.value(x, y+1, z)==1) current.push_back(Pt(x, y+1, z));
if (zsize>z+1 && mask.value(x, y, z+1)==1) current.push_back(Pt(x, y, z+1));
}
// cout<<"make_mask_from_mesh ends"<<endl;
return mask;
}
开发者ID:dungttbk,项目名称:sips,代码行数:50,代码来源:bet2.cpp
示例5: calculateMeans
// calculates the mean for each roi
void RoiMeanCalculation::calculateMeans(volume<float> &actualvolume)
{
if (means.size())
{
for (int i = 0; i < means.size(); i++) means[i] = 0;
for (int z = 0; z < reference.zsize(); z++)
for (int y = 0; y < reference.ysize(); y++)
for (int x = 0; x < reference.xsize(); x++)
{
float value = reference.value(x, y, z);
if (value != 0)
{
int idx = mapping[value];
if (idx >= -1)
means[idx] += actualvolume.value(x, y, z);
}
}
for (int i = 0; i < means.size(); i++)
{
if (counts[i]) means[i] = means[i] / (float)counts[i];
}
}
}
开发者ID:valentinalorenzetti,项目名称:FriendENGINE,代码行数:26,代码来源:masks.cpp
示例6: t1only_special_extract
//extracts profiles in the area around the eyes
vector<double> t1only_special_extract(const volume<float> & t1, const Pt & point, const Vec & n) {
vector<double> resul;
resul.clear();
bool output = true;
const double INNER_DEPTH = 3;
const double OUTER_DEPTH = 100;
Profile pt1;
for (double d = -INNER_DEPTH; d < OUTER_DEPTH; d+=.5)
{
Pt c = point + d * n;
double tmp1 = t1.interpolate(c.X, c.Y, c.Z);
pt1.add (d, tmp1);
}
pt1.init_roi();
//outer skin
double outskin = pt1.last_point_over(pt1.end(), .2);
double check = pt1.last_point_over(outskin - 1.5, .2);
if (outskin - check > 2) output = false;
pt1.set_rroi(outskin);
double inskull = pt1.next_point_under(-INNER_DEPTH, .25);
if (inskull > 5) inskull = 0;
double outskull = pt1.next_point_over(inskull, .35);
resul.push_back(inskull);
resul.push_back(outskull);
resul.push_back(outskin);
return resul;
}
开发者ID:,项目名称:,代码行数:35,代码来源:
示例7: _volume2Sample
void _volume2Sample(svm_model *model, volume<float> &volSample, volume<float> &mask, int sampleSize, float minValue, svm_node * &sample)
{
sample=(struct svm_node *) malloc((sampleSize+1)*sizeof(struct svm_node));
int i=0;
for(int z=0; z < volSample.zsize();z++)
for(int y=0; y < volSample.ysize();y++)
for(int x=0; x < volSample.xsize();x++)
if (mask.value(x,y,z) > minValue)
{
sample[i].index = (i+1);
sample[i].value = volSample.value(x,y,z);
i++;
}
sample[i].index = -1;
}
开发者ID:InstitutoDOr,项目名称:FriendENGINE,代码行数:16,代码来源:svmfuncs.cpp
示例8: draw_segment
void draw_segment(volume<short>& image, const Pt& p1, const Pt& p2)
{
double xdim = (double) image.xdim();
double ydim = (double) image.ydim();
double zdim = (double) image.zdim();
double mininc = min(xdim,min(ydim,zdim)) * .5;
Vec n = p1 - p2;
double d = n.norm();
n.normalize();
for (double i=0; i<=d; i+=mininc)
{
Pt p = p2 + i* n;
image((int) floor((p.X)/xdim +.5),(int) floor((p.Y)/ydim +.5),(int) floor((p.Z)/zdim +.5)) = 1;
}
}
开发者ID:,项目名称:,代码行数:17,代码来源:
示例9: draw_mesh
void draw_mesh(volume<short>& image, const Mesh &m)
{
double xdim = (double) image.xdim();
double ydim = (double) image.ydim();
double zdim = (double) image.zdim();
double mininc = min(xdim,min(ydim,zdim)) * .5;
for (list<Triangle*>::const_iterator i = m._triangles.begin(); i!=m._triangles.end(); i++)
{
Vec n = (*(*i)->get_vertice(0) - *(*i)->get_vertice(1));
double d = n.norm();
n.normalize();
for (double j=0; j<=d; j+=mininc)
{
Pt p = (*i)->get_vertice(1)->get_coord() + j* n;
draw_segment(image, p, (*i)->get_vertice(2)->get_coord());
}
}
}
开发者ID:,项目名称:,代码行数:21,代码来源:
示例10: check_integral
void check_integral(fields &f,
linear_integrand_data &d, const volume &v,
component cgrid)
{
double x1 = v.in_direction_min(d.dx);
double x2 = v.in_direction_max(d.dx);
double y1 = v.in_direction_min(d.dy);
double y2 = v.in_direction_max(d.dy);
double z1 = v.in_direction_min(d.dz);
double z2 = v.in_direction_max(d.dz);
master_printf("Check %d-dim. %s integral in %s cell with %s integrand...",
(x2 - x1 > 0) + (y2 - y1 > 0) + (z2 - z1 > 0),
component_name(cgrid),
v.dim == D3 ? "3d" : (v.dim == D2 ? "2d" :
(v.dim == Dcyl ? "cylindrical"
: "1d")),
(d.c == 1.0 && !d.axy && !d.ax && !d.ay && !d.az
&& !d.axy && !d.ayz && !d.axz) ? "unit" : "linear");
if (0)
master_printf("\n... grid_volume (%g,%g,%g) at (%g,%g,%g) with integral (%g, %g,%g,%g, %g,%g,%g, %g)...\n",
x2 - x1, y2 - y1, z2 - z1,
(x1+x2)/2, (y1+y2)/2, (z1+z2)/2,
d.c, d.ax,d.ay,d.az, d.axy,d.ayz,d.axz, d.axyz);
double sum = real(f.integrate(0, 0, linear_integrand, (void *) &d, v));
if (fabs(sum - correct_integral(v, d)) > 1e-9 * fabs(sum))
abort("FAILED: %0.16g instead of %0.16g\n",
(double) sum, correct_integral(v, d));
master_printf("...PASSED.\n");
}
开发者ID:LeiDai,项目名称:meep,代码行数:31,代码来源:integrate.cpp
示例11: draw_mesh_bis
volume<float> draw_mesh_bis(const volume<float>& image, const Mesh &m)
{
double xdim = (double) image.xdim();
double ydim = (double) image.ydim();
double zdim = (double) image.zdim();
double mininc = min(xdim,min(ydim,zdim)) * .5;
volume<float> res = image;
double max = image.max();
for (list<Triangle*>::const_iterator i = m._triangles.begin(); i!=m._triangles.end(); i++)
{
Vec n = (*(*i)->get_vertice(0) - *(*i)->get_vertice(1));
double d = n.norm();
n.normalize();
for (double j=0; j<=d; j+=mininc)
{
Pt p = (*i)->get_vertice(1)->get_coord() + j* n;
draw_segment_bis(res, p, (*i)->get_vertice(2)->get_coord(), max);
}
}
return res;
}
开发者ID:dungttbk,项目名称:sips,代码行数:22,代码来源:bet2.cpp
示例12: array2Volume
// transforms an array in a volume
void array2Volume(const char *maskFile, float minValue, vector <double> &weightVector, volume<float> &weightVolume)
{
volume<float> mask;
if (maskFile != NULL)
{
string Maskfile = maskFile;
read_volume(mask, Maskfile);
}
weightVolume.reinitialize(mask.xsize(), mask.ysize(), mask.zsize(), 0, true);
weightVolume.copyproperties(mask);
int i = 0;
for(int z=0;z < mask.zsize();z++)
for(int y=0;y < mask.ysize();y++)
for(int x=0;x < mask.xsize();x++)
if (mask.value(x,y,z) > minValue)
{
weightVolume.value(x,y,z) = (float) weightVector[i];
i++;
}
else weightVolume.value(x,y,z) = (float) 0.0;
}
开发者ID:InstitutoDOr,项目名称:FriendENGINE,代码行数:23,代码来源:svmfuncs.cpp
示例13:
void vega::data::hexagonal_prismatic_lattice::fill_volume_cell( volume& v, const prismatic_hexagon_node & node )
{
static const int hex_pixel_offset[8][2] = { {1, 0}, {2, 0}, {0, 1}, {1, 1}, {2, 1}, {3, 1}, {1, 2}, {2, 2}};
math::vector3d vertex = node.get_vertex();
for(size_t h=0;h<8;++h)
{
math::vector3d pixel;
int xx = (int)vertex.x + hex_pixel_offset[h][0];
int yy = (int)vertex.y + hex_pixel_offset[h][1];
int zz = (int)vertex.z;
if( xx >= 0 && yy >= 0 && zz >= 0 && xx < v.get_width() && yy < v.get_height() && zz < v.get_depth() )
{
#ifndef USE_COLOR_MAP
v.set_voxel(xx, yy, zz, (voxel)(node.density * 255.f));
#else
v.set_voxel(xx, yy, zz, (voxel)(node.color.A));
#endif
}
}
}
开发者ID:mihaipopescu,项目名称:VEGA,代码行数:22,代码来源:hexagonal_prismatic_lattice.cpp
示例14: BasisfieldException
void basisfield::AsVolume(volume<float>& vol, FieldIndex fi)
{
if (int(FieldSz_x()) != vol.xsize() || int(FieldSz_y()) != vol.ysize() || int(FieldSz_z()) != vol.zsize()) {
throw BasisfieldException("basisfield::AsVolume:: Matrix size mismatch beween field and volume");
}
if (Vxs_x() != vol.xdim() || Vxs_y() != vol.ydim() || Vxs_z() != vol.zdim()) {
throw BasisfieldException("basisfield::AsVolume:: Voxel size mismatch beween field and volume");
}
if (!coef) {throw BasisfieldException("basisfield::AsVolume: Coefficients undefined");}
if (!UpToDate(fi)) {Update(fi);}
const boost::shared_ptr<NEWMAT::ColumnVector> tmptr = Get(fi);
int vindx=0;
for (unsigned int k=0; k<FieldSz_z(); k++) {
for (unsigned int j=0; j<FieldSz_y(); j++) {
for (unsigned int i=0; i<FieldSz_x(); i++) {
vol(i,j,k) = tmptr->element(vindx++);
}
}
}
}
开发者ID:valentinalorenzetti,项目名称:FriendENGINE,代码行数:22,代码来源:basisfield.cpp
示例15: regionBestVoxels
// returns in `output` the best voxels of `region`
void RegionExtraction::regionBestVoxels(RoiMeanCalculation &reference, volume<float>&values, volume<float>&output, int region, int regionSize, float percentage)
{
vector<roiPoint> roi;
roi.resize(regionSize);
greaterRoiPoint greaterFirst;
int index=0;
// Filter the voxels of the region `region`
for(int z=0;z<values.zsize();z++)
for(int y=0;y<values.ysize();y++)
for(int x=0;x<values.xsize();x++)
{
// get the voxel intensity in reference
int voxelRegion = (int) reference.voxelValues(x,y,z);
// if is the chosen region, records the voxel values (T values or other voxel value)
if (region == voxelRegion)
{
roi[index].value = values.value(x,y,z);
roi[index].roiValue = region;
roi[index].x=x;
roi[index].y=y;
roi[index].z=z;
index++;
}
else values.value(x,y,z)=(float)0.0;
}
// sorts the vector of values in descending order
std::sort(roi.begin(), roi.end(), greaterFirst);
// calculates the cut Index. Remember the voxels descending order
int cutIndex = (int) (roi.size() * percentage + 0.5);
// recording the result in `output`
for (int j=0; j<=cutIndex;j++)
{
output.value(roi[j].x, roi[j].y, roi[j].z) = roi[j].roiValue;
}
}
开发者ID:valentinalorenzetti,项目名称:FriendENGINE,代码行数:42,代码来源:masks.cpp
示例16: correct_integral
static double correct_integral(const volume &v,
const linear_integrand_data &data)
{
direction x = data.dx, y = data.dy, z = data.dz;
double x1 = v.in_direction_min(x);
double x2 = v.in_direction_max(x);
double y1 = v.in_direction_min(y);
double y2 = v.in_direction_max(y);
double z1 = v.in_direction_min(z);
double z2 = v.in_direction_max(z);
return (data.c * integral1(x1,x2,x) * integral1(y1,y2,y) * integral1(z1,z2,z)
+ data.ax * integralx(x1,x2,x) * integral1(y1,y2,y) * integral1(z1,z2,z)
+ data.ay * integral1(x1,x2,x) * integralx(y1,y2,y) * integral1(z1,z2,z)
+ data.az * integral1(x1,x2,x) * integral1(y1,y2,y) * integralx(z1,z2,z)
+ data.axy * integralx(x1,x2,x) * integralx(y1,y2,y) * integral1(z1,z2,z)
+ data.ayz * integral1(x1,x2,x) * integralx(y1,y2,y) * integralx(z1,z2,z)
+ data.axz * integralx(x1,x2,x) * integral1(y1,y2,y) * integralx(z1,z2,z)
+ data.axyz * integralx(x1,x2,x) * integralx(y1,y2,y) * integralx(z1,z2,z)
);
}
开发者ID:LeiDai,项目名称:meep,代码行数:21,代码来源:integrate.cpp
示例17: t1only_co_ext
vector<double> t1only_co_ext(const volume<float> & t1, const Pt & point, const Vec & n) {
vector<double> resul;
resul.clear();
bool output = true;
bool alloutput = true;
const double INNER_DEPTH = 3;
const double OUTER_DEPTH = 60;
Profile pt1;
for (double d = -INNER_DEPTH; d < OUTER_DEPTH; d+=.5)
{
Pt c = point + d * n;
double tmp1 = t1.interpolate(c.X, c.Y, c.Z);
pt1.add (d, tmp1);
}
pt1.init_roi();
//outer skin
double outskin = pt1.last_point_over(pt1.end(), .2);
double check = pt1.last_point_over(outskin - 1.5, .2);
if (outskin - check > 2) outskin = check;
const double OUTER_SKIN = outskin;
pt1.set_rroi(OUTER_SKIN);
double inskull = pt1.next_point_under(pt1.begin(), .25);
pt1.set_lroi(inskull);
if (alloutput)
{
//outer skull
//starting from the skin
double outskull2 = 0;
outskull2 = pt1.last_point_over(outskin, .75);
outskull2 = pt1.last_point_under(outskull2, .20);
//starting from the brain
double minabs = pt1.next_point_under(pt1.begin(), .30);
if (minabs == -500) output = false;
minabs = Max(minabs, -INNER_DEPTH);
double localminabs = minabs;//0
if (output)
{
bool stop = false;
const double lowthreshold = pt1.threshold(.15);
const double upthreshold = pt1.threshold(.60);
int test = 0;
for (vector<pro_pair>::const_iterator i = pt1.v.begin(); i != pt1.v.end(); i++)
{
if (!stop && (*i).abs>=minabs && i!=pt1.v.end() && i!=pt1.v.begin())
{
if ((*i).val>upthreshold) stop = true; //avoid climbing skin
if ((*i).val>lowthreshold) test++;
if ((*i).val<lowthreshold) test--;
if (test < 0) test=0;
if (test == 12) {stop = true;}
if ((*i).val<lowthreshold)
{
localminabs = (*i).abs;
}
}
}
}
double outskull = pt1.next_point_over(localminabs, .15);//.20
if (outskull2 - outskull < -2 || outskull2 - outskull >= 2) {output = false;}
if (outskin - outskull2 < 2) output = false;
if (output)
{
resul.push_back(inskull);
resul.push_back(outskull2);
resul.push_back(outskin);
}
else resul.push_back(outskin);
}
return resul;
}
开发者ID:,项目名称:,代码行数:89,代码来源:
示例18: t1only_write_ext_skull
//writes externall skull computed from image on output.
void t1only_write_ext_skull(volume<float> & output_inskull, volume<float> & output_outskull, volume<float> & output_outskin, const volume<float> & t1, const Mesh & m, const trMatrix & M) {
int glob_counter = 0;
int rem_counter = 0;
const double xdim = t1.xdim();
const double ydim = t1.ydim();
const double zdim = t1.zdim();
double imax = t1.max();
if (imax == 0) imax = 1;
volume<short> meshimage;
copyconvert(t1, meshimage);
meshimage = 0;
draw_mesh(meshimage, m);
for (vector<Mpoint*>::const_iterator i = m._points.begin(); i != m._points.end(); i++)
{
(*i)->data.clear();
double max_neighbour = 0;
const Vec normal = (*i)->local_normal();
const Vec n = Vec(normal.X/xdim, normal.Y/ydim, normal.Z/zdim);
for (list<Mpoint*>::const_iterator nei = (*i)->_neighbours.begin(); nei != (*i)->_neighbours.end(); nei++)
max_neighbour = Max(((**i) - (**nei)).norm(), max_neighbour);
max_neighbour = ceil((max_neighbour)/2);
const Pt mpoint((*i)->get_coord().X/xdim,(*i)->get_coord().Y/ydim,(*i)->get_coord().Z/zdim);
for (int ck = (int)floor(mpoint.Z - max_neighbour/zdim); ck <= (int)floor(mpoint.Z + max_neighbour/zdim); ck++)
for (int cj = (int)floor(mpoint.Y - max_neighbour/ydim); cj <= (int)floor(mpoint.Y + max_neighbour/ydim); cj++)
for (int ci = (int)floor(mpoint.X - max_neighbour/xdim); ci <= (int)floor(mpoint.X + max_neighbour/xdim); ci++)
{
bool compute = false;
const Pt point(ci, cj, ck);
const Pt realpoint(ci*xdim, cj*ydim, ck*zdim);
if (meshimage(ci, cj, ck) == 1)
{
double mindist = 10000;
for (list<Mpoint*>::const_iterator nei = (*i)->_neighbours.begin(); nei != (*i)->_neighbours.end(); nei++)
mindist = Min(((realpoint) - (**nei)).norm(), mindist);
if (mindist >= ((realpoint) - (**i)).norm()) compute = true;
}
if (compute)
{
glob_counter ++;
vector<double> val;
if (!special_case(realpoint, normal, M))
val = t1only_co_ext(t1, point, n);
else
{
val = t1only_special_extract(t1, point, n);
}
if (val.size() == 3)
{
Pt opoint(point.X, point.Y, point.Z);
Vec on(n.X, n.Y, n.Z);
Pt c0 = opoint + val[0]*on;
Pt c1 = opoint + val[1]*on;
Pt c2 = opoint + val[2]*on;
output_inskull((int)floor(c0.X + .5) + infxm,(int) floor(c0.Y + .5) + infym,(int) floor(c0.Z + .5) + infzm) +=1;
output_outskull((int)floor(c1.X + .5) + infxm,(int) floor(c1.Y + .5) + infym,(int) floor(c1.Z + .5) + infzm)+=1;
output_outskin((int)floor(c2.X + .5) + infxm,(int) floor(c2.Y + .5) + infym,(int) floor(c2.Z + .5) + infzm) +=1;
}
else {
rem_counter++;
if (val.size()==1)
{
Pt opoint(point.X, point.Y, point.Z);
Vec on(n.X, n.Y, n.Z);
Pt c0 = opoint + val[0]*on;
output_outskin((int)floor(c0.X + .5) + infxm,(int) floor(c0.Y + .5) + infym,(int) floor(c0.Z + .5) + infzm) +=1;
}
}
}
}
}
if (verbose.value())
{
cout<<" nb of profiles : "<<glob_counter<<endl;
cout<<" removed profiles : "<<100. * rem_counter/(double) glob_counter<<"%"<<endl;
}
}
开发者ID:,项目名称:,代码行数:91,代码来源:
示例19: standard_step_of_computation
double standard_step_of_computation(const volume<float> & image, Mesh & m, const int iteration_number, const double E,const double F, const float addsmooth, const float speed, const int nb_iter, const int id, const int od, const bool vol, const volume<short> & mask){
double xdim = image.xdim();
double ydim = image.ydim();
double zdim = image.zdim();
if (nb_iter % 50 == 0)
{
double l2 = 0;
int counter = 0;
for (vector<Mpoint*>::iterator i = m._points.begin(); i!=m._points.end(); i++ )
{
counter++;
l2 += (*i)->medium_distance_of_neighbours();
}
l = l2/counter;
}
if (nb_iter % 100 == 0)
{
for (vector<Mpoint*>::iterator i = m._points.begin(); i!=m._points.end(); i++)
{
Vec n = (*i)->local_normal();
Pt point = (*i)->get_coord();
Pt ipoint(point.X/xdim, point.Y/ydim, point.Z/zdim);
Vec in(n.X/xdim, n.Y/ydim, n.Z/zdim);
double max = 0;
Pt c_m1 = ipoint + (-1) * in;
double current = image.interpolate((c_m1.X),(c_m1.Y),(c_m1.Z));
for (double i2 = 1; i2 < 150; i2+=2)
{
if (max > .1) break;
Pt c_p = ipoint + i2 * in;
double tmpp = image.interpolate((c_p.X),(c_p.Y),(c_p.Z));
double tmp = (tmpp - current) * 100;
max = Max(max, tmp);
current = tmpp;
if (tmpp > .1) {max = 1; break;}
}
if (max < .1)
{
//There is a problem here for precision mode, since with the copy, no guarantee that data is non-zero size
//even if mesh.cpp operator = is modified to copy data, after retesselate "new" points will have zero size data member
if ( (*i)->data.size() ) (*i)->data.pop_back();
(*i)->data.push_back(1);
}
else
{
if ( (*i)->data.size() ) (*i)->data.pop_back();
(*i)->data.push_back(0);
}
}
}
for (vector<Mpoint*>::iterator i = m._points.begin(); i!=m._points.end(); i++)
{
Vec sn, st, u1, u2, u3, u;
double f2, f3=0;
Vec n = (*i)->local_normal();
Vec dv = (*i)->difference_vector();
double tmp = dv|n;
sn = n * tmp;
st = dv - sn;
u1 = st*.5;
double rinv = (2 * fabs(sn|n))/(l*l);
f2 = (1+tanh(F*(rinv - E)))*0.5;
u2 = f2 * sn * addsmooth;
if ((*i)->data.back() == 0)
{
//main term of skull_extraction
{
Pt point = (*i)->get_coord();
Pt ipoint(point.X/xdim, point.Y/ydim, point.Z/zdim);
Vec in(n.X/xdim, n.Y/ydim, n.Z/zdim);
Pt c_m = ipoint + (-1.) * in;
Pt c_p = ipoint + 1. * in;
double tmp = image.interpolate((c_p.X ),( c_p.Y),(c_p.Z));
double gradient = tmp - image.interpolate((c_m.X),(c_m.Y), (c_m.Z));
double tmp2 = gradient*100;
f3 = max(-1., min(tmp2, 1.));
if (tmp2 >= 0 && tmp2 < .1 && tmp < .1 ) f3 = speed;
if (vol)
{
double tmpvol = mask.interpolate((ipoint.X ),(ipoint.Y),(ipoint.Z));
if (tmpvol > .0)
{
f3 = Max(Max(tmpvol*.5, .1), f3);
f2 = 0;
}
}
//.........这里部分代码省略.........
开发者ID:,项目名称:,代码行数:101,代码来源:
示例20: adjust_initial_mesh
bet_parameters adjust_initial_mesh(const volume<float> & image, Mesh& m, const double & rad = 0., const double xpara=0., const double ypara=0., const double zpara=0.)
{
bet_parameters bp;
double xdim = image.xdim();
double ydim = image.ydim();
double zdim = image.zdim();
double t2, t98, t;
//computing t2 && t98
// cout<<"computing robust min && max begins"<<endl;
bp.min = image.min();
bp.max = image.max();
t2 = image.robustmin();
t98 = image.robustmax();
//t2=32.;
//t98=16121.;
// cout<<"computing robust min && max ends"<<endl;
t = t2 + .1*(t98 - t2);
bp.t98 = t98;
bp.t2 = t2;
bp.t = t;
// cout<<"t2 "<<t2<<" t98 "<<t98<<" t "<<t<<endl;
// cout<<"computing center && radius begins"<<endl;
//finds the COG
Pt center(0, 0, 0);
double counter = 0;
if (xpara == 0. & ypara==0. & zpara==0.)
{
double tmp = t - t2;
for (int k=0; k<image.zsize(); k++)
for (int j=0; j<image.ysize(); j++)
for (int i=0; i<image.xsize(); i++)
{
double c = image(i, j, k ) - t2;
if (c > tmp)
{
c = min(c, t98 - t2);
counter+=c;
center += Pt(c*i*xdim, c*j*ydim, c*k*zdim);
}
}
center=Pt(center.X/counter, center.Y/counter, center.Z/counter);
//cout<<counter<<endl;
// cout<<"cog "<<center.X<<" "<<center.Y<<" "<<center.Z<<endl;
}
else center=Pt(xpara, ypara, zpara);
bp.cog = center;
if (rad == 0.)
{
double radius=0;
counter=0;
double scale=xdim*ydim*zdim;
for (int k=0; k<image.zsize(); k++)
for (int j=0; j<image.ysize(); j++)
for (int i=0; i<image.xsize(); i++)
{
double c = image(i, j, k);
if (c > t)
{
counter+=1;
}
}
radius = pow (.75 * counter*scale/M_PI, 1.0/3.0);
// cout<<radius<<endl;
bp.radius = radius;
}
else (bp.radius = rad);
m.translation(center.X, center.Y, center.Z);
m.rescale (bp.radius/2, center);
// cout<<"computing center && radius ends"<<endl;
//computing tm
// cout<<"computing tm begins"<<endl;
vector<double> vm;
for (int k=0; k<image.zsize(); k++)
for (int j=0; j<image.ysize(); j++)
for (int i=0; i<image.xsize(); i++)
{
double d = image.value(i, j, k);
Pt p(i*xdim, j*ydim, k*zdim);
if (d > t2 && d < t98 && ((p - center)|(p - center)) < bp.radius * bp.radius)
vm.push_back(d);
}
int med = (int) floor(vm.size()/2.);
// cout<<"before sort"<<endl;
nth_element(vm.begin(), vm.begin() + med - 1, vm.end());
//partial_sort(vm.begin(), vm.begin() + med + 1, vm.end());
//double tm = vm[med];
double tm=(*max_element(vm.begin(), vm.begin() + med));
//.........这里部分代码省略.........
开发者ID:dungttbk,项目名称:sips,代码行数:101,代码来源:bet2.cpp
注:本文中的volume类示例由纯净天空整理自Github/MSDocs等源码及文档管理平台,相关代码片段筛选自各路编程大神贡献的开源项目,源码版权归原作者所有,传播和使用请参考对应项目的License;未经允许,请勿转载。 |
请发表评论