void 
Glx::Sphere::sph_interp(Glx::Vector& v1, Glx::Vector& v2, 
			float t, Glx::Vector& res)
{
  double cur[16], inv[16];
  float alpha = acos( dot(v1, v2)/(v1.magnitude()*v2.magnitude()) );
  float delta = alpha*t;
  Glx::Vector i,j,k,tmp, _v1, _v2;

  i.set(v1[0],v1[1],v1[2]);
  i.normalize();

  /* calc the vector ortho to these two, the 'K' vector */
  k = cross(i, v2);

  /* calc a new 'J' vector bases using I and K*/
  j = cross(k,i);
  j.normalize();

  buildMat(cur,i,j,k);
  
  _v1[0] = projection(v1, i);
  _v1[1] = projection(v1, j);
  _v1[2] = projection(v1, k);

  if( Glx::inv4x4(cur, inv) ){
    rot_z(_v1, delta, _v2 );
    Glx::xformVec(_v2, cur, res);
  } else {
    std::cerr << "! INV failed!" << std::endl;
  }
}

Mat nnor::autoCrop(Mat src, int threshold)
{
	Mat verticalProjection = projection(src, VERTICAL);
	Mat horizontalProjection = projection(src, HORIZONTAL);

	int x0, y0, x1, y1;
	for (x0 = 0; x0 < src.cols; x0++)
		if (horizontalProjection.at<int>(x0) > threshold)
			break;
	for (y0 = 0; y0 < src.rows; y0++)
		if (verticalProjection.at<int>(y0) > threshold)
			break;
	for (x1 = src.cols - 1; x1 >= 0; x1--)
		if (horizontalProjection.at<int>(x1) > threshold)
			break;
	for (y1 = src.rows - 1; y1 >= 0; y1--)
		if (verticalProjection.at<int>(y1) > threshold)
			break;

	if ((x1 - x0 <= 0) || (y1 - y0 <= 0))
		return Mat();

	Mat result = src(Rect(x0, y0, x1 - x0 + 1, y1 - y0 + 1));
	return result;
}

//float angle, float ratio, float near, float far, Vector3 &camPos, Vector3 &lookAt, Vector3 &up
void Frustum::Calculate(float depth, const Matrix4 &viewMatrix, const Matrix4 &projectionMatrix)
{
	float zMinimum, r;
	Matrix4 matrix, projection (projectionMatrix), view (viewMatrix);

	
	// Calculate the minimum Z distance in the frustum.
	zMinimum		= -projectionMatrix(4,3) / projectionMatrix(3,3);
	r				= depth / (depth - zMinimum);
	projection(3,3) = r;
	projection(4,3) = -r * zMinimum;

	matrix = view * projection;
	
	// Clear the frustrum & add the 6 planes
	mFrustum.clear();

	Plane near;
	near.a = matrix(1,4) + matrix(1,3);
	near.b = matrix(2,4) + matrix(2,3);
	near.c = matrix(3,4) + matrix(3,3);
	near.d = matrix(4,4) + matrix(4,3);
	mFrustum.push_back(near);

	Plane far;
	far.a = matrix(1,4) - matrix(1,3); 
	far.b = matrix(2,4) - matrix(2,3);
	far.c = matrix(3,4) - matrix(3,3);
	far.d = matrix(4,4) - matrix(4,3);
	mFrustum.push_back(far);

	Plane left;
	left.a = matrix(1,4) + matrix(1,1); 
	left.b = matrix(2,4) + matrix(2,1);
	left.c = matrix(3,4) + matrix(3,1);
	left.d = matrix(4,4) + matrix(4,1);
	mFrustum.push_back(left);

	Plane right;
	right.a = matrix(1,4) - matrix(1,1); 
	right.b = matrix(2,4) - matrix(2,1);
	right.c = matrix(3,4) - matrix(3,1);
	right.d = matrix(4,4) - matrix(4,1);
	mFrustum.push_back(right);

	Plane top;
	top.a = matrix(1,4) - matrix(1,2); 
	top.b = matrix(2,4) - matrix(2,2);
	top.c = matrix(3,4) - matrix(3,2);
	top.d = matrix(4,4) - matrix(4,2);
	mFrustum.push_back(top);

	Plane bottom;
	bottom.a = matrix(1,4) + matrix(1,2);
	bottom.b = matrix(2,4) + matrix(2,2);
	bottom.c = matrix(3,4) + matrix(3,2);
	bottom.d = matrix(4,4) + matrix(4,2);
	mFrustum.push_back(bottom);
		
}

void LeapToCameraCalibrator::resetProjector()
{
	//position = ofVec3f(1.0f,1.0f,1.0f);//cam.getPosition();
	throwRatio = 1.62f;
	lensOffset = ofVec2f(0.0f,0.5f);
    
	ofQuaternion rotation;
	auto rotationQuat = ofQuaternion(rotationX, ofVec3f(1, 0, 0), rotationZ, ofVec3f(0, 0, 1), rotationY, ofVec3f(0, 1, 0));
	ofMatrix4x4 pose = ofMatrix4x4(rotationQuat);
	pose(3,0) = translationX;
	pose(3,1) = translationY;
	pose(3,2) = translationZ;
	projector.setView(pose);
    
	ofMatrix4x4 projection;
	projection(0,0) = throwRatioX;
	projection(1,1) = -throwRatioY;
	projection(2,3) = 1.0f;
	projection(3,3) = 0.0f;
	projection.postMultTranslate(-lensOffsetX, -lensOffsetY, 0.0f);
	projector.setProjection(projection);
    
	projector.setWidth(resolution.x);
	projector.setHeight(resolution.y);
}

// initialize the file for output
void X3DDisplayDevice::write_header(void) {
  fprintf(outfile, "<?xml version=\"1.0\" encoding=\"UTF-8\"?>\n");
  fprintf(outfile, "<!DOCTYPE X3D PUBLIC \"ISO//Web3D//DTD X3D 3.0//EN\"\n");
  fprintf(outfile, "  \"http://www.web3d.org/specifications/x3d-3.0.dtd\">\n");
  fprintf(outfile, "\n");

  // check for special features that require newer X3D versions
  // At present the features that require the latest X3D spec include:
  //  - orthographic camera
  //  - clipping planes
  //  - various special shaders 
  if (projection() == PERSPECTIVE) {
    // if we use a perspective camera, we are compatible with X3D 3.1 spec
    fprintf(outfile, "<X3D version='3.1' profile='Interchange'>\n");
  } else {
    // if we use an orthographic camera, need the newest X3D 3.2 spec
    fprintf(outfile, "<X3D version='3.2' profile='Interchange'>\n");
  }

  fprintf(outfile, "<head>\n");
  fprintf(outfile, "  <meta name='description' content='VMD Molecular Graphics'/>\n");
  fprintf(outfile, "</head>\n");
  fprintf(outfile, "<Scene>\n");
  fprintf(outfile, "<!-- Created with VMD: -->\n");
  fprintf(outfile, "<!-- http://www.ks.uiuc.edu/Research/vmd/ -->\n");

  // export camera definition
  if (projection() == PERSPECTIVE) {
    float vfov = float(2.0*atan2((double) 0.5*vSize, (double) eyePos[2]-zDist));
    if (vfov > VMD_PI)
      vfov=float(VMD_PI); // X3D spec disallows FOV over 180 degrees

    fprintf(outfile, "<Viewpoint description=\"VMD Perspective View\" fieldOfView=\"%g\" orientation=\"0 0 -1 0\" position=\"%g %g %g\" centerOfRotation=\"0 0 0\" />\n",
            vfov, eyePos[0], eyePos[1], eyePos[2]);
  } else {
    fprintf(outfile, "<OrthoViewpoint description=\"VMD Orthographic View\" fieldOfView=\"%g %g %g %g\" orientation=\"0 0 -1 0\" position=\"%g %g %g\" centerOfRotation=\"0 0 0\" />\n",
            -Aspect*vSize/4, -vSize/4, Aspect*vSize/4, vSize/4,
            eyePos[0], eyePos[1], eyePos[2]);
  }

  if (backgroundmode == 1) {
    // emit background sky color gradient
    fprintf(outfile, "<Background skyColor='%g %g %g, %g %g %g, %g %g %g' ", 
            backgradienttopcolor[0], // top pole
            backgradienttopcolor[1],
            backgradienttopcolor[2], 
            (backgradienttopcolor[0]+backgradientbotcolor[0])/2.0f, // horizon
            (backgradientbotcolor[1]+backgradienttopcolor[1])/2.0f, 
            (backgradienttopcolor[2]+backgradientbotcolor[2])/2.0f,
            backgradientbotcolor[0], // bottom pole
            backgradientbotcolor[1],
            backgradientbotcolor[2]);
    fprintf(outfile, "skyAngle='1.5, 3.0' />");
  } else {
    // otherwise emit constant color background sky
    fprintf(outfile, "<Background skyColor='%g %g %g'/>",
            backColor[0], backColor[1], backColor[2]);
  }
  fprintf(outfile, "\n");
}

void RasterImageLayer::loadFile()
{
    GDALDataset *ds = static_cast<GDALDataset*>(GDALOpen(m_filename.toLocal8Bit(), GA_ReadOnly));
    if (ds == nullptr)
    {
        qWarning() << "Error opening file:" << m_filename;
    }
    else
    {
        projection().setGeogCS(new OGRSpatialReference(ds->GetProjectionRef()));
        projection().setProjCS(new OGRSpatialReference(ds->GetProjectionRef()));
        projection().setDomain({-180., -90., 360., 180.});

        std::vector<double> geoTransform(6);
        int xsize = ds->GetRasterXSize();
        int ysize = ds->GetRasterYSize();
        ds->GetGeoTransform(geoTransform.data());
        vertData.resize(4);
        vertData[0] = QVector2D(geoTransform[0], geoTransform[3]);
        vertData[1] = QVector2D(geoTransform[0] + geoTransform[2] * ysize,
                geoTransform[3] + geoTransform[5] * ysize);
        vertData[2] = QVector2D(geoTransform[0] + geoTransform[1] * xsize + geoTransform[2] * ysize,
                geoTransform[3] + geoTransform[4] * xsize + geoTransform[5] * ysize);
        vertData[3] = QVector2D(geoTransform[0] + geoTransform[1] * xsize,
                geoTransform[3] + geoTransform[4] * xsize);
        texData = {{0., 0.}, {0., 1.}, {1., 1.}, {1., 0.}};

        int numBands = ds->GetRasterCount();
        qDebug() << "Bands:" << numBands;

        imData = QImage(xsize, ysize, QImage::QImage::Format_RGBA8888);
        imData.fill(QColor(255, 255, 255, 255));
        // Bands start at 1
        for (int i = 1; i <= numBands; ++i)
        {
            GDALRasterBand *band = ds->GetRasterBand(i);
            switch(band->GetColorInterpretation())
            {
            case GCI_RedBand:
                band->RasterIO(GF_Read, 0, 0, xsize, ysize, imData.bits(),
                               xsize, ysize, GDT_Byte, 4, 0);
                break;
            case GCI_GreenBand:
                band->RasterIO(GF_Read, 0, 0, xsize, ysize, imData.bits() + 1,
                               xsize, ysize, GDT_Byte, 4, 0);
                break;
            case GCI_BlueBand:
                band->RasterIO(GF_Read, 0, 0, xsize, ysize, imData.bits() + 2,
                               xsize, ysize, GDT_Byte, 4, 0);
                break;
            default:
                qWarning() << "Unhandled color interpretation:" << band->GetColorInterpretation();
            }
        }

        GDALClose(ds);
        newFile = true;
    }
}

void FourierProjector::produceSideInfo()
{
    // Zero padding
    MultidimArray<double> Vpadded;
    int paddedDim=(int)(paddingFactor*volumeSize);
    // JMRT: TODO: I think it is a very poor design to modify the volume passed
    // in the construct, it will be padded anyway, so new memory should be allocated
    volume->window(Vpadded,FIRST_XMIPP_INDEX(paddedDim),FIRST_XMIPP_INDEX(paddedDim),FIRST_XMIPP_INDEX(paddedDim),
                   LAST_XMIPP_INDEX(paddedDim),LAST_XMIPP_INDEX(paddedDim),LAST_XMIPP_INDEX(paddedDim));
    volume->clear();
    // Make Fourier transform, shift the volume origin to the volume center and center it
    MultidimArray< std::complex<double> > Vfourier;
    transformer3D.completeFourierTransform(Vpadded,Vfourier);
    ShiftFFT(Vfourier, FIRST_XMIPP_INDEX(XSIZE(Vpadded)), FIRST_XMIPP_INDEX(YSIZE(Vpadded)), FIRST_XMIPP_INDEX(ZSIZE(Vpadded)));
    CenterFFT(Vfourier,true);
    Vfourier.setXmippOrigin();

    // Compensate for the Fourier normalization factor
    double K=(double)(XSIZE(Vpadded)*XSIZE(Vpadded)*XSIZE(Vpadded))/(double)(volumeSize*volumeSize);
    FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(Vfourier)
    DIRECT_MULTIDIM_ELEM(Vfourier,n)*=K;
    Vpadded.clear();
    // Compute Bspline coefficients
    if (BSplineDeg==3)
    {
        MultidimArray< double > VfourierRealAux, VfourierImagAux;
        Complex2RealImag(Vfourier, VfourierRealAux, VfourierImagAux);
        Vfourier.clear();
        produceSplineCoefficients(BSPLINE3,VfourierRealCoefs,VfourierRealAux);
        produceSplineCoefficients(BSPLINE3,VfourierImagCoefs,VfourierImagAux);
        //VfourierRealAux.clear();
        //VfourierImagAux.clear();
    }
    else
        Complex2RealImag(Vfourier, VfourierRealCoefs, VfourierImagCoefs);

    // Allocate memory for the 2D Fourier transform
    projection().initZeros(volumeSize,volumeSize);
    projection().setXmippOrigin();
    transformer2D.FourierTransform(projection(),projectionFourier,false);

    // Calculate phase shift terms
    phaseShiftImgA.initZeros(projectionFourier);
    phaseShiftImgB.initZeros(projectionFourier);
    double shift=-FIRST_XMIPP_INDEX(volumeSize);
    double xxshift = -2 * PI * shift / volumeSize;
    for (size_t i=0; i<YSIZE(projectionFourier); ++i)
    {
        double phasey=(double)(i) * xxshift;
        for (size_t j=0; j<XSIZE(projectionFourier); ++j)
        {
            // Phase shift to move the origin of the image to the corner
            double dotp = (double)(j) * xxshift + phasey;
            sincos(dotp,&DIRECT_A2D_ELEM(phaseShiftImgB,i,j),&DIRECT_A2D_ELEM(phaseShiftImgA,i,j));
        }
    }
}

void
Panner::draw_the_box ( int tx, int ty, int tw, int th )
{
    draw_box();

    Fl_Image *i = 0;

    if ( _bg_image && ( _bg_image->h() != th || projection() != _bg_image_projection ) )
    {
        if ( _bg_image_scaled )
            delete _bg_image;
        else
            ((Fl_Shared_Image*)_bg_image)->release();

        _bg_image = 0;
    }

    if ( ! _bg_image )
    {
        if ( projection() == POLAR )
        {
            if ( th <= 92 )
                i = Fl_Shared_Image::get( PIXMAP_PATH "/non-mixer/panner-sphere-92x92.png" );
            else if ( th <= 502 )
                i = Fl_Shared_Image::get( PIXMAP_PATH "/non-mixer/panner-sphere-502x502.png" );
            else if ( th <= 802 )
                i = Fl_Shared_Image::get( PIXMAP_PATH "/non-mixer/panner-sphere-802x802.png" );
        }
        else
        {
            if ( th <= 92 )
                i = Fl_Shared_Image::get( PIXMAP_PATH "/non-mixer/panner-plane-92x92.png" );
            else if ( th <= 502 )
                i = Fl_Shared_Image::get( PIXMAP_PATH "/non-mixer/panner-plane-502x502.png" );
            else if ( th <= 802 )
                i = Fl_Shared_Image::get( PIXMAP_PATH "/non-mixer/panner-plane-802x802.png" );
        }

        if ( i && i->h() != th )
        {
            Fl_Image *scaled = i->copy( th, th );

            _bg_image = scaled;
            _bg_image_scaled = true;
        }
        else
        {
            _bg_image = i;
            _bg_image_scaled = false;
        }
    }
    
    _bg_image_projection = projection();

    if ( _bg_image )
        _bg_image->draw( tx, ty );
}

int main(int argc, char** argv) {
    if (2>argc) {
        std::cerr << "Usage: " << argv[0] << "obj/model.obj" << std::endl;
        return 1;
    }

    float *zbuffer = new float[width*height];
    for (int i=width*height; i--; zbuffer[i] = -std::numeric_limits<float>::max());

    TGAImage frame(width, height, TGAImage::RGB);
    lookat(eye, center, up);
    viewport(width/8, height/8, width*3/4, height*3/4);
    projection(-1.f/(eye-center).norm());
    light_dir = proj<3>((Projection*ModelView*embed<4>(light_dir, 0.f))).normalize();

    for (int m=1; m<argc; m++) {
        model = new Model(argv[m]);
        Shader shader;
        for (int i=0; i<model->nfaces(); i++) {
            for (int j=0; j<3; j++) {
                shader.vertex(i, j);
            }
            triangle(shader.varying_tri, shader, frame, zbuffer);
        }
        delete model;
    }
    frame.flip_vertically(); // to place the origin in the bottom left corner of the image
    frame.write_tga_file("framebuffer.tga");

    delete [] zbuffer;
    return 0;
}

int main() {

    std::string inputfile = "tests/008-ch.osc.gz";

    options_t options;

    std::shared_ptr<reprojection> projection(reprojection::create_projection(PROJ_SPHERE_MERC));
    options.projection = projection;

    auto out_test = std::make_shared<test_output_t>(options);
    osmdata_t osmdata(std::make_shared<dummy_slim_middle_t>(), out_test, options.projection);

    boost::optional<std::string> bbox;
    parse_osmium_t parser(bbox, true, &osmdata);

    parser.stream_file(inputfile, "");

    assert_equal(out_test->node.added, 0);
    assert_equal(out_test->node.modified, 1176);
    assert_equal(out_test->node.deleted, 16773);
    assert_equal(out_test->way.added, 0);
    assert_equal(out_test->way.modified, 161);
    assert_equal(out_test->way.deleted, 4);
    assert_equal(out_test->rel.added, 0);
    assert_equal(out_test->rel.modified, 11);
    assert_equal(out_test->rel.deleted, 1);

    return 0;
}

/*
	Function for setting up and running the fluid simulation kernels
*/
void solveFluid(struct Configuration* config) { 
	// Jacobi settings
	float alpha = -1.0f;
	// Should use alpha -1 here, but this gives nicer results
	//float alpha = -(1.0f/invhalfgridscale);
	float rbeta = 0.25;
	int iterations = 100;

	// grid scaling. this is currently not used
    if(rank == 0){
        float gridscale = 1.0f;
        float invgridscale = 1.0f/gridscale;
        float invhalfgridscale = 0.5f/gridscale; 

        // Timstep value
        float timestep = 0.05f;

        // Emitter settings
        float amount = 2.0f;
        float radius = 0.5*config->N/10.0f;
        float emitterposx = config->N/2;
        float emitterposy = config->N/3;

        // buoyancy settings
        float bdiry = 1.0f;
        float bdirx = 0.0f;
        float bstrength = 0.1f;

        // advection settings
        float veldamp = 0.01f;	


        // Velocity advection
        float *tmp = config->velx0; config->velx0 = config->velx; config->velx = tmp;
        float *tmp1 = config->vely0; config->vely0 = config->vely; config->vely = tmp1;
        advect(config->N, config->velx, config->velx0, config->velx0, config->vely0, timestep, veldamp, 1);
        advect(config->N, config->vely, config->vely0, config->velx0, config->vely0, timestep, veldamp, 2);

        // Density advection
        float *tmp2 = config->dens0; config->dens0 = config->dens; config->dens = tmp2;
        advect(config->N, config->dens, config->dens0, config->velx, config->vely, timestep, 0.0f, 0);

        // Add density and density buoyancy
        addDensity(config->N, config->dens, timestep, emitterposx, emitterposy, radius, amount);
        addDensityBuoyancy(config->N, config->velx, config->vely, config->dens, bdirx, bdiry, bstrength, timestep);

        // Divergence calculation
        divergence(config->N, config->velx, config->vely, config->div);	

        // Pressure jacobi calculation. First set pres array to zero as initial guess
        setMem(config->N, config->pres);
    }

	jacobi(iterations);
	
    if(rank == 0){
        // Calculate projection
        projection(config->N, config->velx, config->vely, config->pres);	
    }
}

int main(int argc, char *argv[]) {
  char *srcdir = getenv("srcdir");

  if (srcdir == NULL) {
    std::cerr << "$srcdir not set!\n";
    return 1;
  }

  std::string inputfile = std::string(srcdir) + std::string("/tests/test_multipolygon.osm");

  options_t options;
  boost::shared_ptr<reprojection> projection(new reprojection(PROJ_SPHERE_MERC));
  options.projection = projection;

  boost::shared_ptr<test_output_t> out_test(new test_output_t(options));
  osmdata_t osmdata(boost::make_shared<test_middle_t>(), out_test);

  parse_xml2_t parser(0, false, projection, 0, 0, 0, 0);

  int ret = parser.streamFile(inputfile.c_str(), 0, &osmdata);
  if (ret != 0) {
    return ret;
  }

  assert_equal(out_test->sum_ids,       73514L);
  assert_equal(out_test->num_nodes,       353L);
  assert_equal(out_test->num_ways,        140L);
  assert_equal(out_test->num_relations,    40L);
  assert_equal(out_test->num_nds,         495L);
  assert_equal(out_test->num_members,     146L);

  return 0;
}

void PreviewWindow::pickPrimitive()
{
    const std::vector<std::shared_ptr<Primitive>> &prims = _scene->primitives();

    std::shared_ptr<Camera> cam = _scene->camera();

    MatrixStack::set(VIEW_STACK, cam->transform());
    MatrixStack::set(PROJECTION_STACK, projection());

    _fbo->bind();
    _fbo->attachDepthBuffer(*_depthBuffer);
    _fbo->attachTexture(*_screenBuffer, 0);
    _fbo->selectAttachments(1);
    _fbo->setReadBuffer(RT_ATTACHMENT0);

    glViewport(0, 0, width(), height());

    glDisable(GL_MULTISAMPLE);
    glClearColor(1.0f, 1.0f, 1.0f, 1.0f);
    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);

    _solidShader->bind();
    renderMeshes(*_solidShader, [this](uint32 i) {
        _solidShader->uniformF(
            "Color",
            ((i      ) & 0xFF)/255.0f,
            ((i >>  8) & 0xFF)/255.0f,
            ((i >> 16) & 0xFF)/255.0f,
            ((i >> 24) & 0xFF)/255.0f
        );
        return true;
    });

/*!

 */
bool
Segment2D::onSegmentWeakly( const Vector2D & p ) const
{
    Vector2D proj = projection( p );

    return ( proj.isValid()
             && p.equalsWeakly( proj ) );

#if 0
    Vector2D o = origin();
    Vector2D t = terminal();

    const double buf = ( allow_on_terminal
                         ? EPSILON
                         : 0.0 );

    if ( std::fabs( ( t - o ).outerProduct( p - o ) ) < EPSILON )
    {
        if ( std::fabs( o.x - t.x ) < EPSILON )
        {
            return ( ( o.y - buf < p.y && p.y < t.y + buf )
                     || ( t.y - buf < p.y && p.y < o.y + buf ) );
        }
        else
        {
            return ( ( o.x - buf < p.x && p.x < t.x + buf )
                     || ( t.x - buf < p.x && p.x < o.x + buf) );
        }
    }

    return false;
#endif
}

vector<int> Dist::projectionOnestep(vector<double> Y, vector<double> & e, vector< vector<double> > & X)
{
    e.resize(Y.size(), 0);
    vector< vector<double> > Xtmp, Xtmp2;
    Xtmp=X;
    vector<double> beta(Xtmp.size(),0);
    vector<bool> significant(Xtmp.size(),false);

    vector<int> idx, idx2;

    for(int i=0; i<X.size(); i++) idx.push_back(i);

    projectionSolu(Y,e,beta,significant,Xtmp);

    Xtmp2.clear();
    idx2.clear();
    for(int i=0; i<significant.size(); i++) {
        if(significant[i]) {
            Xtmp2.push_back(Xtmp[i]);
            idx2.push_back(idx[i]);
        }
    }

    if(Xtmp2.size()==0) {
        e=Y;
    } else {
        projection(Y,e,Xtmp2);
    }


    return idx2;
}

void updateProj(AttributePtr<Attribute_Camera> camera)
{
	DirectX::XMFLOAT4X4& projection = *(DirectX::XMFLOAT4X4*)&camera->mat_projection;
	float& zFar			= camera->zFar;
	float& zNear		= camera->zNear;
	float& aspectRatio	= camera->aspectRatio;
	float& fieldOfView	= camera->fieldOfView;

	memset(&projection, 0, sizeof(projection));

	projection(0, 0) = 1/(aspectRatio*(tan(fieldOfView/2)));
	projection(1, 1) = 1/(tan(fieldOfView/2));
	projection(2, 2) = zFar/(zFar - zNear);
	projection(2, 3) = 1.0f;
	projection(3, 2) = (-zNear * zFar)/(zFar - zNear);
}

    Shader::Shader(const View& view, SDL_Window& window)
        : mOriginalViewport()
        , mProgram(loadProgram("assets/shaders/basic.glvs", "assets/shaders/basic.glfs"))
        , mProjectionLocation(glGetUniformLocation(mProgram, "te_ProjectionMatrix"))
        , mModelViewLocation(glGetUniformLocation(mProgram, "te_ModelViewMatrix"))
        , mProjection()
    {
        glGetIntegerv(GL_VIEWPORT, mOriginalViewport.data());

        int width = 0, height = 0;
        SDL_GetWindowSize(&window, &width, &height);
        FloatRect viewport = view.getViewport();
        glViewport((int)(viewport.x * width),
                   (int)(height - (height * (viewport.h + viewport.y))),
                   (int)(viewport.w * width),
                   (int)(viewport.h * height));

        glUseProgram(mProgram);

        FloatRect lens = view.getLens();
        glm::mat4 projection(glm::ortho<GLfloat>(lens.x, lens.x + lens.w, lens.y + lens.h, lens.y, -Z, Z));
        if (mProjectionLocation == -1) { throw std::runtime_error("te_ProjectionMatrix: not a valid program variable."); }
        glUniformMatrix4fv(mProjectionLocation, 1, GL_FALSE, glm::value_ptr(projection));

        if (mModelViewLocation == -1) { throw std::runtime_error{ "te_ModelViewMatrix: not a valid program variable." }; }
        glUniformMatrix4fv(mModelViewLocation, 1, GL_FALSE, glm::value_ptr(glm::mat4()));
    }

void VeryFastView::set
  (Voxel *v,int d)
  {
    //  if(v->color.a<0.7)
    //    return;
    if(!parent->inIso(v->pos))
      return;



    Pos2D p2=projection(v->pos);
    //  cdebug("size:"<<mViewSpace.size());
    std::map<Pos2D,std::pair<int,Voxel*> >::iterator i=mViewSpace.find(p2);
    //  cdebug("size:"<<mViewSpace.size());

    //  CTRACE;
    //  cdebug(p2<<"//"<<v->pos);

    if(i!=mViewSpace.end())
      {
        if(i->second.first>d)
          {
            i->second=std::make_pair(d,v);
            //         cdebug("set:"<<i->second.first<<":"<<d);
            return;
          }
        //      cdebug(i->second.first<<":"<<d);
        return;
      }
    //  cdebug("set2:"<<d);
    mViewSpace[p2]=std::make_pair(d,v);
  }

int main() {
    srand((unsigned) time(NULL));

    buf = initBuf(BUFSIZE, BLKSIZE);

    initData();

    //输出R的内容
    //printData(R_START);
    //输出S的内容
    //printData(S_START);

    //找到R.A = 40
    int addr_find_R_A = findfirst(R_START, 40);
    printData(addr_find_R_A);

    //找到S.C = 60
    int addr_find_S_C = findfirst(S_START, 60);
    printData(addr_find_S_C);


    //R.A投影
    int projection_R_A = projection(R_START);
    printData(projection_R_A);
    printf("NUM = %d\n", BLKNUM);
    //R.A连接S.C
    int r_join_s = join(R_START, S_START);
    //printData(r_join_s);

    //printf("%d\n", BLKNUM);
    return 0;
}

void TextEngine::draw( const kvs::Vec3& p, const std::string& text, kvs::ScreenBase* screen ) const
{
    GLdouble model[16]; kvs::OpenGL::GetModelViewMatrix( model );
    GLdouble proj[16]; kvs::OpenGL::GetProjectionMatrix( proj );
    GLint view[4]; kvs::OpenGL::GetViewport( view );
    GLdouble winx = 0, winy = 0, winz = 0;
    kvs::OpenGL::Project( p.x(), p.y(), p.z(), model, proj, view, &winx, &winy, &winz );

    kvs::OpenGL::WithPushedAttrib attrib( GL_ALL_ATTRIB_BITS );
    attrib.disable( GL_TEXTURE_1D );
    attrib.disable( GL_TEXTURE_2D );
//    attrib.enable( GL_TEXTURE_2D );
    attrib.disable( GL_TEXTURE_3D );
    attrib.enable( GL_DEPTH_TEST );
    attrib.enable( GL_BLEND );
    {
        kvs::OpenGL::SetBlendFunc( GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA );
        kvs::OpenGL::WithPushedMatrix modelview( GL_MODELVIEW );
        modelview.loadIdentity();
        {
            kvs::OpenGL::WithPushedMatrix projection( GL_PROJECTION );
            projection.loadIdentity();
            {
                const GLint left = view[0];
                const GLint top = view[1];
                const GLint right = view[0] + view[2];
                const GLint bottom = view[1] + view[3];
                kvs::OpenGL::SetOrtho( left, right, bottom, top, 0, 1 );
                kvs::OpenGL::Translate( 0, 0, -winz );
                m_font.draw( kvs::Vec2( winx, top - ( winy - bottom ) ), text );
            }
        }
    }
}