C++ Vec::dot方法代码示例

//alg for calculating shading. Calls diffuseShade and SpecularShade
Color RT_lights(Figure* obj, const Ray& ray, const Vec& thePoint, const Vec& normal)
{
    Color c = Color();
    for (list<Light*>::iterator iterator = lightList.begin(), end = lightList.end(); iterator != end; ++iterator) {
        Light* theLight = *iterator;

        pair<double, Figure*> inter = nearestIntersection(Ray(Vec(thePoint), theLight->getPos()), MIN_T / SHAD_RES, 1.0, EPSILON, false);
        if (inter.first <= 0) {
            Vec* toLight = thePoint.normalize(theLight->getPos());
            double dotProduct = toLight->dot(normal);
            Vec* subt = (thePoint.sub(theLight->getPos()));
            double dist = abs(subt->getMag());
            Color dif = diffuseShade(obj, theLight, max(dotProduct, 0.0), dist);
            Vec* Q = normal.scale(normal.dot(*toLight));
            Vec* S = Q->sub(*toLight);
            Vec* S2 = S->scale(2.0);
            Vec* R = toLight->add(*S2);
            Vec* norm = ray.getNorm();
            Vec* scaledNorm = norm->scale(-1.0);
            dotProduct = max(0.0, R->dot(*scaledNorm));
            Color spec = specularShade(obj, normal, theLight, *toLight, pow(dotProduct, obj->getShininess()), ray, dist);
            c = c.add(dif.add(spec));
            delete(toLight);
            delete(Q);
            delete(S);
            delete(R);
            delete(S2);
            delete(subt);
            delete(norm);
            delete(scaledNorm);
        }
    }
    return c;
}

inline double bssrdf(const Vec& xi, const Vec& ni, const Vec& wi, const Vec& xo, const Vec& no, const Vec& wo, const int j) {
	// distance
	const Vec xoxi = xo - xi;
	const double r = xoxi.len();

	// modified normal
	const Vec ni_s = (xoxi.normalized()) % ((ni % xoxi).normalized());

	// directions of ray sources
	const double nnt = 1.0 / eta, ddn = -wi.dot(ni);
	const Vec wr = (wi * -nnt - ni * (ddn * nnt + sqrt(1.0 - nnt * nnt * (1.0 - ddn * ddn)))).normalized();
	const Vec wv = wr - ni_s * (2.0 * wr.dot(ni_s));

	// distance to real sources
	const double cos_beta = -sqrt((r * r - xoxi.dot(wr) * xoxi.dot(wr)) / (r * r + de[j] * de[j]));
	double dr;
	const double mu0 = -no.dot(wr);
	if (mu0 > 0.0) {
		dr = sqrt((D[j] * mu0) * ((D[j] * mu0) - de[j] * cos_beta * 2.0) + r * r);
	} else {
		dr = sqrt(1.0 / (3.0 * sigma_t[j] * 3.0 * sigma_t[j]) + r * r);
	}

	// distance to virtual source
	const Vec xoxv = xo - (xi + ni_s * (2.0 * A * de[j]));
	const double dv = xoxv.len();

	// BSSRDF
	const double result = Sp_d(xoxi, wr, dr, no, j) - Sp_d(xoxv, wv, dv, no, j);

	// clamping to zero
	return (result < 0.0) ? 0.0 : result;
}

 double intersect(const Ray &r) const { // returns distance, 0 if nohit
     Vec op = p - r.o; // Solve t^2*d.d + 2*t*(o-p).d + (o-p).(o-p)-R^2 = 0
     double t, eps = 1e-4, b = op.dot(r.d),
               det = b * b - op.dot(op) + rad * rad;
     if (det < 0)
         return 0;
     else
         det = sqrt(det);
     return (t = b - det) > eps ? t : ((t = b + det) > eps ? t : 0);
 }

	inline double intersect(const Ray& r) const {
		// ray-sphere intersection returns the distance
		const Vec op = p - r.o;
		double t, b = op.dot(r.d), det = b * b - op.dot(op) + rad * rad;
		if (det < 0) {
			return 1e20;
		}
		else {
			det = sqrt(det);
		}
		return (t = b - det) > 1e-4 ? t : ((t = b + det) > 1e-4 ? t : 1e20);
	}

// directional dipole
// --------------------------------
inline double Sp_d(const Vec& x, const Vec& w, const double& r, const Vec& n, const int j) {
	// evaluate the profile
	const double s_tr_r = sigma_tr[j] * r;
	const double s_tr_r_one = 1.0 + s_tr_r;
	const double x_dot_w = x.dot(w);
	const double r_sqr = r * r;

	const double t0 = Cp_norm * (1.0 / (4.0 * PI * PI)) * exp(-s_tr_r) / (r * r_sqr);
	const double t1 = r_sqr / D[j] + 3.0 * s_tr_r_one * x_dot_w;
	const double t2 = 3.0 * D[j] * s_tr_r_one * w.dot(n);
	const double t3 = (s_tr_r_one + 3.0 * D[j] * (3.0 * s_tr_r_one + s_tr_r * s_tr_r) / r_sqr * x_dot_w) * x.dot(n);

	return t0 * (Cp * t1 - Ce * (t2 - t3));
}

Number Sphere::do_intersect(const Vec & start_p, const Vec & dir) const {
    Vec p = start_p - center;
    Vec x = (-dir.dot(p)) * dir;
    Vec y = p + x;
    auto y_norm = y.norm();
    if(y_norm > radius)
        return -1;
    auto x_dot = x.dot(dir);

    auto s = std::sqrt(radius * radius - y_norm * y_norm);
    auto t_norm = x_dot - s;
    if(t_norm < 0)
        t_norm += 2 * s;
    return t_norm; // may < 0, which means no intersection
}

Vec radiance(const Ray &r, int depth, unsigned short *Xi,int E=1){
  double t;                               // distance to intersection
  int id=0;                               // id of intersected object
  if (!intersect(r, t, id)) return Vec(); // if miss, return black
  const Sphere &obj = spheres[id];        // the hit object
  Vec x=r.o+r.d*t, n=(x-obj.p).norm(), nl=n.dot(r.d)<0?n:n*-1, f=obj.c;
  double p = f.x>f.y && f.x>f.z ? f.x : f.y>f.z ? f.y : f.z; // max refl
  if (++depth>5||!p) if (erand48(Xi)<p) f=f*(1/p); else return obj.e*E;
  if (obj.refl == DIFF){                  // Ideal DIFFUSE reflection
    double r1=2*M_PI*erand48(Xi), r2=erand48(Xi), r2s=sqrt(r2);
    Vec w=nl, u=((fabs(w.x)>.1?Vec(0,1):Vec(1))%w).norm(), v=w%u;
    Vec d = (u*cos(r1)*r2s + v*sin(r1)*r2s + w*sqrt(1-r2)).norm();

    // Loop over any lights
    Vec e;
    for (int i=0; i<numSpheres; i++){
      const Sphere &s = spheres[i];
      if (s.e.x<=0 && s.e.y<=0 && s.e.z<=0) continue; // skip non-lights
      
      Vec sw=s.p-x, su=((fabs(sw.x)>.1?Vec(0,1):Vec(1))%sw).norm(), sv=sw%su;
      double cos_a_max = sqrt(1-s.rad*s.rad/(x-s.p).dot(x-s.p));
      double eps1 = erand48(Xi), eps2 = erand48(Xi);
      double cos_a = 1-eps1+eps1*cos_a_max;
      double sin_a = sqrt(1-cos_a*cos_a);
      double phi = 2*M_PI*eps2;
      Vec l = su*cos(phi)*sin_a + sv*sin(phi)*sin_a + sw*cos_a;
      l.norm();
      if (intersect(Ray(x,l), t, id) && id==i){  // shadow ray
        double omega = 2*M_PI*(1-cos_a_max);
        e = e + f.mult(s.e*l.dot(nl)*omega)*M_1_PI;  // 1/pi for brdf
      }
    }
    
    return obj.e*E+e+f.mult(radiance(Ray(x,d),depth,Xi,0));
  } else if (obj.refl == SPEC)              // Ideal SPECULAR reflection
    return obj.e + f.mult(radiance(Ray(x,r.d-n*2*n.dot(r.d)),depth,Xi));
  Ray reflRay(x, r.d-n*2*n.dot(r.d));     // Ideal dielectric REFRACTION
  bool into = n.dot(nl)>0;                // Ray from outside going in?
  double nc=1, nt=1.5, nnt=into?nc/nt:nt/nc, ddn=r.d.dot(nl), cos2t;
  if ((cos2t=1-nnt*nnt*(1-ddn*ddn))<0)    // Total internal reflection
    return obj.e + f.mult(radiance(reflRay,depth,Xi));
  Vec tdir = (r.d*nnt - n*((into?1:-1)*(ddn*nnt+sqrt(cos2t)))).norm();
  double a=nt-nc, b=nt+nc, R0=a*a/(b*b), c = 1-(into?-ddn:tdir.dot(n));
  double Re=R0+(1-R0)*c*c*c*c*c,Tr=1-Re,P=.25+.5*Re,RP=Re/P,TP=Tr/(1-P);
  return obj.e + f.mult(depth>2 ? (erand48(Xi)<P ?   // Russian roulette
    radiance(reflRay,depth,Xi)*RP:radiance(Ray(x,tdir),depth,Xi)*TP) :
    radiance(reflRay,depth,Xi)*Re+radiance(Ray(x,tdir),depth,Xi)*Tr);
}

Vec receivedRadiance(const Ray &r, int depth, bool flag) {
    double t;                                   // Distance to intersection
    int id = 0;                                 // id of intersected sphere

    if ( !intersect(r, t, id) ) return Vec();   // if miss, return black
    const Sphere &obj = spheres[id];            // the hit object

    Vec x = r.o + r.d*t;                        // The intersection point
    Vec o = (Vec() - r.d).normalize();          // The outgoing direction (= -r.d)

    Vec n = (x - obj.p).normalize();            // The normal direction
    if ( n.dot(o) < 0 ) n = n*-1.0;

    /*
    Tips

    1. Other useful quantities/variables:
    Vec Le = obj.e;                             // Emitted radiance
    const BRDF &brdf = obj.brdf;                // Surface BRDF at x

    2. Call brdf.sample() to sample an incoming direction and continue the recursion
    */

    Vec Le = obj.e;                             // Emitted radiance
    const BRDF &brdf = obj.brdf;                // Surface BRDF at x

    const int rrDepth = 5;
    const double survivalProbability = 0.9;

    double p = 1.0;

    if (depth > rrDepth)
        p = survivalProbability;

    if (rng() < p)
    {
        Vec o_i;
        double pdf;
        brdf.sample(n, o, o_i, pdf);
        Ray y(x, o_i.normalize());
        Vec reflected = receivedRadiance(y, depth + 1, false)
            .mult(brdf.eval(n, o, o_i))
            * (n.dot(o_i)
                / (pdf*p));
        return Le + reflected;
    } 
    return Le;
}

Vec radiance(const Ray &r, int depth, unsigned short *Xi) {
	double t;                               // distance to intersection
	int id = 0;                               // id of intersected object
	if(!intersect(r, t, id)) return Vec(); // if miss, return black
	const Sphere &obj = spheres[id];        // the hit object
	Vec x = r.o + r.d*t, n = (x - obj.p).norm(), nl = n.dot(r.d) < 0 ? n : n*-1, f = obj.c;
	double p = f.x > f.y && f.x>f.z ? f.x : f.y > f.z ? f.y : f.z; // max refl
	if(depth > 255) return obj.e;
	if(++depth > 5) if(erand48(Xi) < p) f = f*(1 / p); else return obj.e; //R.R.
	if(obj.refl == DIFF) {                  // Ideal DIFFUSE reflection
		double r1 = 2 * M_PI*erand48(Xi), r2 = erand48(Xi), r2s = sqrt(r2);
		Vec w = nl, u = ((fabs(w.x) > .1 ? Vec(0, 1) : Vec(1)) % w).norm(), v = w%u;
		Vec d = (u*cos(r1)*r2s + v*sin(r1)*r2s + w*sqrt(1 - r2)).norm();
		return obj.e + f.mult(radiance(Ray(x, d), depth, Xi));
	}
	else if(obj.refl == SPEC)            // Ideal SPECULAR reflection
		return obj.e + f.mult(radiance(Ray(x, r.d - n * 2 * n.dot(r.d)), depth, Xi));
	Ray reflRay(x, r.d - n * 2 * n.dot(r.d));     // Ideal dielectric REFRACTION
	bool into = n.dot(nl) > 0;                // Ray from outside going in?
	double nc = 1, nt = 1.5, nnt = into ? nc / nt : nt / nc, ddn = r.d.dot(nl), cos2t;
	if((cos2t = 1 - nnt*nnt*(1 - ddn*ddn)) < 0)    // Total internal reflection
		return obj.e + f.mult(radiance(reflRay, depth, Xi));
	Vec tdir = (r.d*nnt - n*((into ? 1 : -1)*(ddn*nnt + sqrt(cos2t)))).norm();
	double a = nt - nc, b = nt + nc, R0 = a*a / (b*b), c = 1 - (into ? -ddn : tdir.dot(n));
	double Re = R0 + (1 - R0)*c*c*c*c*c, Tr = 1 - Re, P = .25 + .5*Re, RP = Re / P, TP = Tr / (1 - P);
	return obj.e + f.mult(depth > 2 ? (erand48(Xi) < P ?   // Russian roulette
		radiance(reflRay, depth, Xi)*RP : radiance(Ray(x, tdir), depth, Xi)*TP) :
		radiance(reflRay, depth, Xi)*Re + radiance(Ray(x, tdir), depth, Xi)*Tr);
}

int main() {
    double y_vals[] = {-1.5, 2, -2.5};
    double z_vals[] = {3, -2, 1};
    Vec<double> zeroes(3);                             // Vec size 3 (entries initialize to zero)
    Vec<double> x = Vec<double>::constantVec(3, 2.5);  // Vec size 3 with all entries set to 2.5
    Vec<double> y = Vec<double>(y_vals, 3);
    Vec<double> z(3);
    z.setEntries(z_vals, 3);
    Vec<int> ix(x);

    cout << "zeroes = " << zeroes << endl;
    cout << "x = " << x << endl;
    cout << "y = " << y << endl;
    cout << "z = " << z << endl;
    cout << "ix = " << ix << endl;
    cout << "z[0] = " << z[0] << ", z[1] = " << z[1] << ", z[2] = " << z[2] << endl;
    cout << "3.5 * x = " << (3.5 * x) << endl;
    cout << "x / 3.5 = " << (x / 3.5) << endl;
    cout << "x + y = " << (x + y) << endl;
    cout << "x - y = " << (x - y) << endl;
    cout << "x.concatenate(y) = " << x.concatenate(y) << endl;
    cout << "x.dot(y) = " << x.dot(y) << endl;
    cout << "x.cross(y) = " << x.cross(y) << endl;
    cout << "x.norm() = " << x.norm() << endl;
    cout << "x.unit_vector() = " << x.unit_vector() << endl;
    cout << "ix.norm() = " << ix.norm() << endl;
    cout << "ix.norm<double>() = " << ix.norm<double>() << endl;
    cout << "ix.unit_vector<double>() = " << ix.unit_vector<double>() << endl;
    cout << "scalar_triple_product(x, y, z) = " 
         << Vec<double>::scalar_triple_product(x, y, z) << endl;
    cout << "vector_triple_product(x, y, z) = " 
         << Vec<double>::vector_triple_product(x, y, z) << endl;
}

bool UVSphere::intersect(const Ray& r, float tmin, float tmax, float time, HitRecord& record)const{
  Vec temp = r.origin() - center;
  float a = r.direction().dot(r.direction());
  float b = 2 * r.direction().dot(temp);
  float c = temp.dot(temp) - radius * radius;

  float discriminant = b*b - 4*a*c;
  if(discriminant > 0){
    discriminant = sqrt(discriminant);
    float t = (-b -discriminant) /(2*a);
    if(t < tmin) t = (-b + discriminant) / (2*a);
    if(t < tmin || t > tmax) return false;
    
    record.t = t;
    record.hit_p = r.origin() + r.direction() * t;
    record.reflect = reflectionCoef;
    record.transparency = refractionCoef;
    Vec n = record.normal = normalize(r.origin() + r.direction()*t - center);
    
    //calculate UV coordinates
    float theta = acos(n.y);
    float phi = atan2(n.z, n.x);
    if(phi < 0) phi += M_PI * 2;
    record.uv = Vec(phi/(2* M_PI), (M_PI - theta)/M_PI, 0);
    record.hit_tex = tex;
    return true;
  }
  return false;
}

  //----------------------------------------------------------------------
  void BLCSSS::rwm_draw_chunk(int chunk){
    clock_t start = clock();
    const Selector &inc(m_->coef().inc());
    int nvars = inc.nvars();
    Vec full_nonzero_beta = m_->beta();   // only nonzero components
    // Compute information matrix for proposal distribution.  For
    // efficiency, also compute the log-posterior of the current beta.
    Vec mu(inc.select(pri_->mu()));
    Spd siginv(inc.select(pri_->siginv()));
    double original_logpost = dmvn(full_nonzero_beta, mu, siginv, 0, true);

    const std::vector<Ptr<BinomialRegressionData> > &data(m_->dat());
    int nobs = data.size();

    int full_chunk_size = compute_chunk_size();
    int chunk_start = chunk * full_chunk_size;
    int elements_remaining = nvars - chunk_start;
    int this_chunk_size = std::min(elements_remaining, full_chunk_size);
    Selector chunk_selector(nvars, false);
    for(int i = chunk_start; i< chunk_start + this_chunk_size; ++i) {
      chunk_selector.add(i);
    }

    Spd proposal_ivar = chunk_selector.select(siginv);

    for(int i = 0; i < nobs; ++i){
      Vec x = inc.select(data[i]->x());
      double eta = x.dot(full_nonzero_beta);
      double prob = plogis(eta);
      double weight = prob * (1-prob);
      VectorView x_chunk(x, chunk_start, this_chunk_size);
      // Only upper triangle is accessed.  Need to reflect at end of loop.
      proposal_ivar.add_outer(x_chunk, weight, false);
      int yi = data[i]->y();
      int ni = data[i]->n();
      original_logpost += dbinom(yi, ni, prob, true);
    }
    proposal_ivar.reflect();
    VectorView beta_chunk(full_nonzero_beta, chunk_start, this_chunk_size);
    if(tdf_ > 0){
      beta_chunk = rmvt_ivar_mt(
          rng(), beta_chunk, proposal_ivar / rwm_variance_scale_factor_, tdf_);
    }else{
      beta_chunk = rmvn_ivar_mt(
          rng(), beta_chunk, proposal_ivar / rwm_variance_scale_factor_);
    }

    double logpost = dmvn(full_nonzero_beta, mu, siginv, 0, true);
    Vec full_beta(inc.expand(full_nonzero_beta));
    logpost += m_->log_likelihood(full_beta, 0, 0, false);
    double log_alpha = logpost - original_logpost;
    double logu = log(runif_mt(rng()));
    ++rwm_chunk_attempts_;
    if(logu < log_alpha){
      m_->set_beta(full_nonzero_beta);
      ++rwm_chunk_successes_;
    }
    clock_t end = clock();
    rwm_chunk_times_ += double(end - start) / CLOCKS_PER_SEC;
  }

/** Cast reflected ray
 ** @param ray incoming ray
 ** @param pt intersection point
 ** @param normal normal of pt
 ** @return reflected ray
 **/
Ray reflect(Ray ray, Vec pt, Vec normal){
  Vec v = ray.getOrig() - pt; v = v.normalize();
  Vec dir = normal * (2 * v.dot(normal)) - v; dir = dir.normalize(); // reflected direction
  Ray t(pt, dir, 0, 9999, ReflectedRay);
  if (debugMode) {
    printf("reflected ray: origin = "); t.getOrig().print();
    printf("reflected ray: direction = "); t.getDir().print();
  }
  return t;
}

void HardSpherePressureTracker::update_collision(Box &box, AtomID a1, AtomID a2,
        flt time, Vec delta_p) {
    if (isnan(t0)) {
        t0 = time;
        lastt = t0;
        return;
    }
    Vec dr = box.diff(a1->x, a2->x);
    collisionsum += dr.dot(delta_p);
    Ksum += kinetic_energy_com(*atoms);
    Ncollisions++;
    lastt = time;
};

/// Kinetic energy relative to center of mass
// K = ½Σmᵢvᵢ² - ½(Σmᵢvᵢ)² / (Σmᵢ)
flt kinetic_energy_com(AtomGroup &atoms) {
    flt Ksum = 0;
    Vec mvsum = Vec::Zero();
    flt msum = 0;
    for (uint i = 0; i < atoms.size(); i++) {
        Atom &a = atoms[i];
        Vec mv = a.m * a.v;
        Ksum += mv.dot(a.v);
        mvsum += mv;
        msum += a.m;
    }
    return (Ksum - mvsum.squaredNorm() / msum) / 2;
};