C++ ScalarField类代码示例

Volume* createFloatFieldVolumeFromVolume(Volume *base_volume)
{
    int width  = base_volume->width();
    int height = base_volume->height();
    int depth  = base_volume->depth();

    std::vector<AbstractScalarField*> fields;

    for(int m=0; m < base_volume->numberOfMaterials(); m++)
    {
        float *data = new float[width*height*depth];

        for(int d=0; d < depth; d++){
            for(int h=0; h < height; h++){
                for(int w=0; w < width; w++){
                    data[d*width*height + h*width + w] = (float)base_volume->valueAt(w+0.5,h+0.5,d+0.5, m);
                }
            }
        }

        ScalarField<float> *floatField = new ScalarField<float>(data,width,height,depth);
        floatField->setName(base_volume->getMaterial(m)->name() + "_asFloat");
        fields.push_back(floatField);
    }

    cleaver::Volume *floatVolume = new Volume(fields);
    floatVolume->setName(base_volume->name());

    return floatVolume;
}

TEST_F(FieldFETest, scalar) {
	FE_P_disc<1,1,3> fe1;
	FE_P_disc<1,2,3> fe2;
	FE_P_disc<1,3,3> fe3;
    ScalarField field;

    dh->distribute_dofs(fe1, fe2, fe3);
    field.set_fe_data(dh, &map1, &map2, &map3, &v);
    field.set_time(0.0);

    dof_values[0] = 1;
    dof_values[1] = 2;
    dof_values[2] = 3;

    vector<double> values(3);

    // test values at vertices of the triangle
    field.value_list( { { 1, 1, 5 }, { 4, 0, 5 }, { 2, 3, 5 } }, mesh->element_accessor(0), values );
    EXPECT_DOUBLE_EQ( dof_values[0], values[0] );
    EXPECT_DOUBLE_EQ( dof_values[1], values[1] );
    EXPECT_DOUBLE_EQ( dof_values[2], values[2] );

    // test value at barycenter
    EXPECT_DOUBLE_EQ( (dof_values[0]+dof_values[1]+dof_values[2])/3, field.value({ 7./3, 4./3, 5 }, mesh->element_accessor(0)) );
}

ScalarField<float>* createFloatFieldFromScalarField(AbstractScalarField *scalarField)
{
    int w = (int)(scalarField->bounds().size.x);
    int h = (int)(scalarField->bounds().size.y);
    int d = (int)(scalarField->bounds().size.z);

    int wh = w*h;
    int whd = wh*d;

    float *data = new float[whd];

    for(int k=0; k < d; k++)
    {
        for(int j=0; j < h; j++)
        {
            for(int i=0; i < w; i++)
            {
                data[k*wh + j*w + i] = (float)scalarField->valueAt(i+0.5, j+0.5, k+0.5);
            }
        }
    }

    ScalarField<float> *floatField = new ScalarField<float>(data, w, h, d);
    floatField->setBounds(scalarField->bounds());

    return floatField;
}

ScalarField<double>* createDoubleFieldFromScalarField(AbstractScalarField *scalarField)
{
    int w = (int)(scalarField->bounds().size.x);
    int h = (int)(scalarField->bounds().size.y);
    int d = (int)(scalarField->bounds().size.z);

    int wh = w*h;
    int whd = wh*d;

    double *data = new double[whd];

    for(int k=0; k < d; k++)
    {
        for(int j=0; j < h; j++)
        {
            for(int i=0; i < w; i++)
            {
                data[k*wh + j*w + i] = scalarField->valueAt(i+0.5, j+0.5, k+0.5);
            }
        }
    }

    ScalarField<double> *doubleField = new ScalarField<double>(data, w, h, d);
    doubleField->setBounds(scalarField->bounds());

    return doubleField;
}

bool ChunkedPointCloud::enableScalarField()
{
    ScalarField* currentInScalarField = getCurrentInScalarField();

	if (!currentInScalarField)
	{
		//if we get there, it means that either the caller has forgot to create
		//(and assign) a scalar field to the cloud, or that we are in a compatibility
		//mode with old/basic behaviour: a unique SF for everything (input/output)

        //we look for any already existing "default" scalar field 
		m_currentInScalarFieldIndex = getScalarFieldIndexByName("Default");
		if (m_currentInScalarFieldIndex < 0)
		{
            //if not, we create it
            m_currentInScalarFieldIndex = addScalarField("Default");
            if (m_currentInScalarFieldIndex < 0) //Something went wrong
                return false;
		}

        currentInScalarField = getCurrentInScalarField();
        assert(currentInScalarField);
	}

	//if there's no output scalar field either, we set this new scalar field as output also
	if (!getCurrentOutScalarField())
		m_currentOutScalarFieldIndex = m_currentInScalarFieldIndex;

	return currentInScalarField->resize(m_points->capacity());
}

int ChunkedPointCloud::addScalarField(const char* uniqueName, bool isStrictlyPositive)
{
    //we don't accept two SF with the same name!
    if (getScalarFieldIndexByName(uniqueName)>=0)
        return -1;

	//create requested scalar field
    ScalarField* sf = new ScalarField(uniqueName,isStrictlyPositive);
	if (!sf)
		return -1;

	//we don't want 'm_scalarFields' to grow by 50% each time! (default behavior of std::vector::push_back)
	try
	{
		m_scalarFields.push_back(sf);
	}
	catch (.../*const std::bad_alloc&*/) //out of memory
	{
		sf->release();
		return -1;
	}

	sf->link();

	return (int)m_scalarFields.size()-1;
}

int ChunkedPointCloud::addScalarField(const char* uniqueName)
{
    //we don't accept two SF with the same name!
    if (getScalarFieldIndexByName(uniqueName) >= 0)
        return -1;

	//create requested scalar field
    ScalarField* sf = new ScalarField(uniqueName);
	if (!sf || (size() && !sf->resize(size())))
	{
		//Not enough memory!
		if (sf)
			sf->release();
		return -1;
	}

	try
	{
		//we don't want 'm_scalarFields' to grow by 50% each time! (default behavior of std::vector::push_back)
		m_scalarFields.resize(m_scalarFields.size()+1);
	}
	catch (std::bad_alloc) //out of memory
	{
		sf->release();
		return -1;
	}
	
	m_scalarFields.back() = sf;
	sf->link();

	return (int)m_scalarFields.size()-1;
}

void deriv_yy(const ScalarField<T, S>& f, ScalarField<S, S>& out) {
  const GridParams& p = f.grid_params();
  assert(p == out.grid_params());
  assert(p.ny >= 3);

  const double d = 1./(p.eps_y*p.eps_y);

  for (int i = 0; i < p.nx; ++i) {
    for (int j = 0; j < p.ny; ++j) {
      T a, b, c;
      if (j == 0) {
        a = f(i,j);
        b = f(i,j+1);
        c = f(i,j+2);
      } else if (j == p.ny-1) {
        a = f(i,j-2);
        b = f(i,j-1);
        c = f(i,j);
      } else {
        a = f(i,j-1);
        b = f(i,j);
        c = f(i,j+1);
      }
      out(i,j) = (a - 2.*b + c) * d;
    }
  }
}

void deriv_xx(const ScalarField<T, S>& f, ScalarField<S, S>& out) {
  const GridParams& p = f.grid_params();
  assert(p == out.grid_params());
  assert(p.nx >= 3);

  const double d = 1./(p.eps_x*p.eps_x);

  for (int i = 0; i < p.nx; ++i) {
    for (int j = 0; j < p.ny; ++j) {
      T a, b, c;
      if (i == 0) {
        a = f(i,j);
        b = f(i+1,j);
        c = f(i+2,j);
      } else if (i == p.nx-1) {
        a = f(i-2,j);
        b = f(i-1,j);
        c = f(i,j);
      } else {
        a = f(i-1,j);
        b = f(i,j);
        c = f(i+1,j);
      }
      out(i,j) = (a - 2.*b + c) * d;
    }
  }
}

ScalarField Real(const complexScalarField& C)
{   const GridInfo& gInfo = C->gInfo;
    ScalarField R;
    nullToZero(R, gInfo);
    callPref(eblas_daxpy)(gInfo.nr, C->scale, (const double*)C->dataPref(false), 2, R->dataPref(false), 1);
    R->scale = 1;
    return R;
}

void apply_flow(const ScalarField<ElemT, ExprT>& phi, const DoubleField& u_x, const DoubleField& u_y, ScalarField<ExprT, ExprT>& out) {
  const GridParams& gp = phi.grid_params();
  assert(u_x.grid_params() == gp && u_y.grid_params() == gp && out.grid_params() == gp);

  for (int i = 0; i < gp.nx; ++i) {
    for (int j = 0; j < gp.ny; ++j) {
      auto xy = gp.to_xy(i, j);
      out(i,j) = phi.eval_xy(xy.first - u_x(i,j), xy.second - u_y(i,j));
    }
  }
}

void deriv_y(const ScalarField<T, S>& f, ScalarField<S, S>& g_y) {
  const GridParams& p = f.grid_params();
  assert(p == g_y.grid_params());

  // derivatives in y direction
  for (int i = 0; i < p.nx; ++i) {
    for (int j = 0; j < p.ny - 1; ++j) {
      g_y(i,j) = (f(i,j+1) - f(i,j)) * (1./p.eps_y);
    }
    // copy over last column
    g_y(i,p.ny-1) = g_y(i,p.ny-2);
  }
}

bool ChunkedPointCloud::renameScalarField(int index, const char* newName)
{
	if (getScalarFieldIndexByName(newName) < 0)
	{
		ScalarField* sf = getScalarField(index);
		if (sf)
		{
			sf->setName(newName);
			return true;
		}
	}
	return false;
}

void deriv_x(const ScalarField<T, S>& f, ScalarField<S, S>& g_x) {
  const GridParams& p = f.grid_params();
  assert(p == g_x.grid_params());

  // derivatives in x direction
  for (int i = 0; i < p.nx - 1; ++i) {
    for (int j = 0; j < p.ny; ++j) {
      g_x(i,j) = (f(i+1,j) - f(i,j)) * (1./p.eps_x);
    }
  }
  // copy over last row
  for (int j = 0; j < p.ny; ++j) {
    g_x(p.nx-1,j) = g_x(p.nx-2,j);
  }
}

int main(int argc, char *argv[])
{

  char str[200];
  int i;


  if(argc<7 )
    {
      std::cout<<"Usage: dtiFA  xsize ysize zsize input_dtfile output_fa_file switch_endian_input(0/1)\n";
      exit(0);
    }
  
 
  int Xdim= atoi(argv[1]);
  int Ydim= atoi(argv[2]);
  int Zdim= atoi(argv[3]);
  
  TensorField dti;
  
  dti.init(Xdim,Ydim,Zdim);
  
  int endian_be=0; endian_be= atoi(argv[6]);

  if(endian_be)
    {
      //std::cout<<"Reading in big endian format\n";
      dti.importRawFieldFromFile_BE(argv[4]);
    }
  else
    {
      dti.importRawFieldFromFile(argv[4]);
    }

  //std::cout<<"DTI read\n";

  ScalarField fa;
  fa.init(Xdim,Ydim,Zdim);
  // fa.setVoxelSize(0,Xres); fa.setVoxelSize(1,Yres); fa.setVoxelSize(2,Zres);
  dti.ComputeFA(&fa);
  //std::cout<<"Smooth FA computed\n";
  
  sprintf(str,"%s",argv[5]);
  fa.exportRawFieldToFile(str);
  
  return 0;

}