projection.cpp

Go to the documentation of this file.

/* projection.cpp :  Stereographic and related projection computations
 * Copyright (C) 2020  Daniel Haley
 * 
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 * 
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 * 
 * You should have received a copy of the GNU General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */

#include "atomprobe/projection/projection.h"

#include "atomprobe/helper/misc.h"

#include "atomprobe/helper/aptAssert.h"

#include <cmath>
#include <sstream>

//GNU Scientific library, for root solving algorithms
#include <gsl/gsl_errno.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_roots.h>

namespace AtomProbe
{

//return residual of tan(theta) - sin(eta)/(m+cos(eta)) =0 
struct SphericProjectionParams
{
    double theta;
    double focusDist;
};

double SphericProjectionEqn(double eta, void *p)
{
    SphericProjectionParams *params=(SphericProjectionParams*)p;

    return tan(params->theta) - sin(eta)/(params->focusDist+cos(eta));
}



bool GnomonicProjection::toPlanar(float theta ,float phi, 
                    float &fx, float &fy) const 
{
    //Only defined for +ve half of sphere
    if(phi >= M_PI/2.0)
        return false;

    //Two equations important here.
    // 1) sphere equation, in polar coordinates
    // 2) plane @ z=1

    //find the point on the sphere surface
    Point3D pSphere;
    pSphere.setISOSpherical(1,theta,phi);

    ASSERT(pSphere[2] >0);

    //Draw line from centre (0,0,0) of sphere to point (P') on plane 
    // which intersects pSphere
    // this implies that [P']=t[pSphere], but we know P'_z = 1
    // so t = 1/pSphere_z
    float t;

    t=1.0/pSphere[2];
    
    Point3D pPlane=pSphere*t;

    fx= pPlane[0];
    fy= pPlane[1];
    
    return true; 
}


bool GnomonicProjection::toAzimuthal(float fx, float fy,
                    float &theta, float &phi) const
{
    //theta is formed by the triangle that extends to the plane, with
    // a radius of r
    // theta is given by atan2
    float r= sqrt(fx*fx + fy*fy);
    
    theta=atan(r);
    phi=atan2f(fy,fx);

    //transform is always valid for any real finite fx,fy pair
    // but is not unique (0,0)
    return true;
}

void GnomonicProjection::scaleDown(float flightLength,
            float detX,float detY,float &scaledX, float &scaledY) const
{
    ASSERT(flightLength > 0);


    //similar triangles. scaledRadius/1 = radius/flightLength
    // 1 due to radius of unit circle
    scaledX= detX/flightLength;
    scaledY= detY/flightLength;
}

void GnomonicProjection::scaleUp(float flightLength,float scaledX,float scaledY,
            float &realX, float &realY) const
{
    ASSERT(flightLength > 0);

    //similar triangles. scaledRadius/2 = radius/flightLength
    // 1 due to radius of unit circle
    realX=scaledX*flightLength;
    realY=scaledY*flightLength;
}

bool StereographicProjection::toPlanar(float theta ,float phi, 
                    float &fx, float &fy) const 
{
    //Multiplied by 2, as this is the triangle length
    // to the rear of the sphere (from front)
    float r = 2.0*cos(theta);
    fx= r*cos(phi);
    fy= r*sin(phi);

    return true;
}

void StereographicProjection::scaleDown(float flightLength,
            float detX,float detY,float &scaledX, float &scaledY) const
{
    ASSERT(flightLength > 0);

    //similar triangles. scaledRadius/2 = radius/flightLength
    // 2 due to diameter of unit circle
    scaledX= 2.0*detX/flightLength;
    scaledY= 2.0*detY/flightLength;
}

void StereographicProjection::scaleUp(float flightLength,float scaledX, float scaledY,
            float &realX,float &realY) const
{
    ASSERT(flightLength > 0);

    //similar triangles. scaledRadius/2 = radius/flightLength
    // 2 due to diameter of unit circle
    realX = scaledX*flightLength/2.0;
    realY = scaledY*flightLength/2.0;
}

bool StereographicProjection::toAzimuthal(float fx, float fy,
                    float &theta, float &phi) const
{
    // phi is given by atan2
    phi=atan2f(fy,fx);
    
    //theta is formed by the triangle that extends to the plane, with
    // a radius of r, base length 2 (2*unit sphere size)
    float r= sqrt(fx*fx + fy*fy);
    theta=atan(r/2.0f);
    
    //transform is always valid for any real finite fx,fy pair
    // but is not unique at (0,0)
    return true;
}


float StereographicProjection::thetaToEta(float theta)  const
{
    //TODO: Fixme - this is a repeat of MFSP code
    //actual equation is tan(theta) = sin(eta)/(1+cos(eta))
    
    //Use brent/brent-dekker solver on 0,pi to solve
    const gsl_root_fsolver_type *T = gsl_root_fsolver_brent;
    gsl_root_fsolver *s = gsl_root_fsolver_alloc (T);

    gsl_function f;
    f.function  =&SphericProjectionEqn;
    SphericProjectionParams m;
    m.focusDist=1;
    m.theta=theta;
    f.params = &m;

    unsigned int iter=0;
    double x_lo=0.0;
    double x_hi=M_PI/2;
    double r; 

    gsl_root_fsolver_set (s, &f, x_lo, x_hi);
    int status;
    do
    {
        iter++;
        gsl_root_fsolver_iterate (s);
        r = gsl_root_fsolver_root (s);
        x_lo = gsl_root_fsolver_x_lower (s);
        x_hi = gsl_root_fsolver_x_upper (s);
        status = gsl_root_test_interval (x_lo, x_hi,
                           0, 0.001);

    } while (status == GSL_CONTINUE && iter < 20);  

    
    gsl_root_fsolver_free(s);   

    if(status != GSL_SUCCESS)
        return -1.0f;
    
    return r;
}

ModifiedFocusSphericProjection::ModifiedFocusSphericProjection(float focus)
{
    setFocus(focus);
}

void ModifiedFocusSphericProjection::setFocus(float focus)
{
    //To be a far side projection, this has to be >1, to lie outside the unit circle,
    // if < 1, it is inside the circle. Cannot be <0, or will not project to plane
//  ASSERT(focus>0);
    focusDist = focus;
}


float ModifiedFocusSphericProjection::etaToTheta(float eta) const
{
    //Eta is the spherical angle, which should be between 0 and PI radiians
    ASSERT(eta <= M_PI && eta>=0);

    //the eqn is tan(theta) = sin(eta)/( cos(eta) + m)
    return atan( sin(eta) / (cos(eta) + focusDist));
}

float ModifiedFocusSphericProjection::thetaToEta(float theta) const
{

    //actual equation is tan(theta) = sin(eta)/(m+cos(eta))
    
    //Use brent/brent-dekker solver on 0,pi to solve
    const gsl_root_fsolver_type *T = gsl_root_fsolver_brent;
    gsl_root_fsolver *s = gsl_root_fsolver_alloc (T);

    gsl_function f;
    f.function  =&SphericProjectionEqn;
    SphericProjectionParams m;
    m.focusDist=focusDist;
    m.theta=theta;
    f.params = &m;

    unsigned int iter=0;
    double x_lo=0.0;
    double x_hi=M_PI/2;
    double r; 

    gsl_root_fsolver_set (s, &f, x_lo, x_hi);
    int status;
    do
    {
        iter++;
        gsl_root_fsolver_iterate (s);
        r = gsl_root_fsolver_root (s);
        x_lo = gsl_root_fsolver_x_lower (s);
        x_hi = gsl_root_fsolver_x_upper (s);
        status = gsl_root_test_interval (x_lo, x_hi,
                           0, 0.001);

    } while (status == GSL_CONTINUE && iter < 20);  

    
    gsl_root_fsolver_free(s);   

    if(status != GSL_SUCCESS)
        return -1.0f;
    
    return r;
}


float ModifiedFocusSphericProjection::getMaxFOV() const
{
    if(focusDist >1.0)
    {
        //Derived from : Line through focus at (-focusDist,0), unit sphere (centered at 0)
        // and projection line at x=1, tangent to sphere
        return acos(-1.0f/focusDist);
    }
    else
    {
        //Focus between 0 and 1
        //This is the angle formed backwards from the orgin to the focus
        return acos(-focusDist);

    }

}

float ModifiedFocusSphericProjection::getFOVRadius(float eta) const
{
    //convert to angle at focus point
    float theta=etaToTheta(eta);

    return (1+focusDist)*tan(theta);
}

bool ModifiedFocusSphericProjection::getProjectionPt(float theta, float &f) const
{
    //If does not intersect sphere, then do not allow
    if(theta >=getMaxFOV())
        return false;

    //intersection of line through focus (-focusDist,0)
    // and projection line at x=1
    f=(1.0f + focusDist)*tan(theta);

    return true;
}


bool ModifiedFocusSphericProjection::toPlanar(float theta, float phi, float &fx, float &fy) const
{
    float fRadius;
    if(!getProjectionPt(theta,fRadius) )
        return false;


    fx = fRadius*cos(phi);
    fy = fRadius*sin(phi);
    
    return true;
}


bool ModifiedFocusSphericProjection::toAzimuthal(float fx, float fy, float &theta, float &phi) const
{
    float radius;

    radius=sqrt(fx*fx + fy*fy);
    theta = atan(radius/(1+focusDist));

    phi=atan2(fy,fx);

    if(theta > getMaxFOV())
        return false;

    return true;
}

void ModifiedFocusSphericProjection::scaleDown(float flightLength,
            float detX,float detY,float &scaledX,float &scaledY) const
{
    ASSERT(flightLength > 0);

    //similar triangles. scaledRadius/2 = radius/flightLength
    // 2 due to diameter of unit circle
    scaledX= (1.0+focusDist)*detX/flightLength;
    scaledY= (1.0+focusDist)*detY/flightLength;
}

void ModifiedFocusSphericProjection::scaleUp(float flightLength,float scaledX,float scaledY,
            float &realX,float &realY) const
{
    ASSERT(flightLength > 0);

    //similar triangles. scaledRadius/2 = radius/flightLength
    // 2 due to diameter of unit circle
    realX= scaledX*flightLength/(1.0+focusDist);
    realY= scaledY*flightLength/(1.0+focusDist);
}


#ifdef DEBUG

bool testAzimuthalProjection();

bool testProjection()
{
    if(!testAzimuthalProjection())
        return false;

    return true;
}

bool testAzimuthalProjection()
{
    {
    GnomonicProjection gp;

    //Test projection onto sphere
    float fx,fy;
    fx=0.0;
    fy=1.0;
    
    float theta,phi;
    gp.toAzimuthal(fx,fy,theta,phi);
    TEST(tolEqual(theta,(float)(45.0f*M_PI/180.0f),0.1f), "toAzimuthal failed");

    //Test projection from sphere to plane
    theta=45.0f*M_PI/180.0f;
    phi=0;

    fx=fy=0;
    gp.toPlanar(theta,phi,fx,fy);
    TEST(tolEqual(fx,1.0f,0.1f),"toPlanar failed (x)");
    TEST(tolEqual(fy,0.0f,0.1f),"toPlanar failed (y)");
    }

    {
    ModifiedFocusSphericProjection mSp1(1.0f);
    TEST(tolEqual(mSp1.etaToTheta(0.875129) ,0.43756f,0.01f),"Projection Eqn");
    }

    {
    ModifiedFocusSphericProjection mSp(2.0f);

    float theta,phi;

    mSp.toAzimuthal(3,0,theta,phi);

    TEST(tolEqual(phi,0.0f,0.1f),"On x axis");
    TEST(tolEqual(theta,(float)(45.0f*M_PI/180.0f),0.1f),"Should form 1-1-sqrt(2) triangle");
    TEST(theta < mSp.getMaxFOV(),"Inside FOV");

    double roundTripTheta;
    roundTripTheta=mSp.etaToTheta(mSp.thetaToEta(0.1));
    TEST(tolEqual(roundTripTheta,0.1, 0.01),"Round-trip coordinate transform for modified stereographic");
    }

    

    return true;
}



#endif
}