coords.h

Go to the documentation of this file.

 
#ifndef OAKFIELD_STIMULUS_COORDS_H
#define OAKFIELD_STIMULUS_COORDS_H
 
#include "oakfield/field.h"
#include "oakfield/field_patch.h"
 
#include <math.h>
#include <stdbool.h>
 
#ifdef __cplusplus
extern "C" {
#endif
 
#define STIMULUS_COORD_EPS 1.0e-12
 
typedef enum SimStimulusCoordMode {
    SIM_STIMULUS_COORD_AXIS = 0,  
    SIM_STIMULUS_COORD_ANGLE,     
    SIM_STIMULUS_COORD_RADIAL,    
    SIM_STIMULUS_COORD_POLAR,     
    SIM_STIMULUS_COORD_AZIMUTH,   
    SIM_STIMULUS_COORD_ELLIPTIC,  
    SIM_STIMULUS_COORD_SEPARABLE, 
    SIM_STIMULUS_COORD_SPIRAL     
} SimStimulusCoordMode;
 
typedef enum SimStimulusCoordAxis {
    SIM_STIMULUS_AXIS_X = 0, 
    SIM_STIMULUS_AXIS_Y = 1  
} SimStimulusCoordAxis;
 
typedef enum SimStimulusSeparableMode {
    SIM_STIMULUS_SEPARABLE_MULTIPLY = 0, 
    SIM_STIMULUS_SEPARABLE_ADD = 1       
} SimStimulusSeparableMode;
 
typedef struct SimStimulusCoordConfig {
    SimStimulusCoordMode mode;        
    SimStimulusCoordAxis axis;        
    SimStimulusSeparableMode combine; 
    double angle;                     
    double origin_x;                  
    double origin_y;                  
    double spacing_x;                 
    double spacing_y;                 
    double center_x;                  
    double center_y;                  
    double velocity_x;                
    double velocity_y;                
    double ellipse_u;                 
    double ellipse_v;                 
    double spiral_arms;               
    double spiral_pitch;              
    double spiral_phase;              
    double spiral_angular_velocity;   
} SimStimulusCoordConfig;
 
typedef struct SimStimulusCoordRow {
    size_t x0;            
    size_t y0;            
    size_t width;         
    double x;             
    double y;             
    double x_step;        
    double y_step;        
    double sample_x;      
    double sample_y;      
    double sample_x_step; 
    double sample_y_step; 
} SimStimulusCoordRow;
 
static inline void sim_stimulus_coord_sample_xy(const SimStimulusCoordConfig *coord, double x,
                                                double y, double t, double *out_x, double *out_y) {
    if (out_x == NULL || out_y == NULL) {
        return;
    }
 
    double sample_x = x;
    double sample_y = y;
    if (coord != NULL) {
        sample_x -= coord->velocity_x * t;
        sample_y -= coord->velocity_y * t;
    }
 
    *out_x = sample_x;
    *out_y = sample_y;
}
 
static inline void sim_stimulus_coord_centered_xy(const SimStimulusCoordConfig *coord, double x,
                                                  double y, double t, double *out_dx,
                                                  double *out_dy) {
    if (out_dx == NULL || out_dy == NULL) {
        return;
    }
 
    double sample_x = x;
    double sample_y = y;
    double cx = 0.0;
    double cy = 0.0;
    if (coord != NULL) {
        sim_stimulus_coord_sample_xy(coord, x, y, t, &sample_x, &sample_y);
        cx = coord->center_x;
        cy = coord->center_y;
    }
 
    *out_dx = sample_x - cx;
    *out_dy = sample_y - cy;
}
 
static inline void sim_stimulus_coord_rotate_xy(double x, double y, double angle, double *out_u,
                                                double *out_v) {
    if (out_u == NULL || out_v == NULL) {
        return;
    }
 
    double s = sin(angle);
    double c = cos(angle);
    *out_u = x * c + y * s;
    *out_v = -x * s + y * c;
}
 
static inline void sim_stimulus_coord_elliptic_local(const SimStimulusCoordConfig *coord, double dx,
                                                     double dy, double *out_u, double *out_v) {
    if (out_u == NULL || out_v == NULL) {
        return;
    }
 
    double ur = dx;
    double vr = dy;
    double au = 1.0;
    double av = 1.0;
 
    if (coord != NULL) {
        sim_stimulus_coord_rotate_xy(dx, dy, coord->angle, &ur, &vr);
        au = coord->ellipse_u;
        av = coord->ellipse_v;
    }
 
    if (fabs(au) <= STIMULUS_COORD_EPS) {
        au = 1.0;
    }
    if (fabs(av) <= STIMULUS_COORD_EPS) {
        av = 1.0;
    }
 
    *out_u = ur / au;
    *out_v = vr / av;
}
 
static inline void sim_stimulus_coord_elliptic_polar(const SimStimulusCoordConfig *coord, double dx,
                                                     double dy, double *out_r, double *out_theta) {
    double u = 0.0;
    double v = 0.0;
 
    sim_stimulus_coord_elliptic_local(coord, dx, dy, &u, &v);
 
    if (out_r != NULL) {
        *out_r = hypot(u, v);
    }
    if (out_theta != NULL) {
        *out_theta = atan2(v, u);
    }
}
 
static inline void sim_stimulus_coord_polar(const SimStimulusCoordConfig *coord, double x, double y,
                                            double t, double *out_r, double *out_theta) {
    double dx = 0.0;
    double dy = 0.0;
 
    sim_stimulus_coord_centered_xy(coord, x, y, t, &dx, &dy);
 
    if (out_r != NULL) {
        *out_r = hypot(dx, dy);
    }
    if (out_theta != NULL) {
        *out_theta = atan2(dy, dx);
    }
}
 
static inline void sim_stimulus_coord_normalize(SimStimulusCoordConfig *coord) {
    if (coord == NULL) {
        return;
    }
    if (coord->mode < SIM_STIMULUS_COORD_AXIS || coord->mode > SIM_STIMULUS_COORD_SPIRAL) {
        coord->mode = SIM_STIMULUS_COORD_AXIS;
    }
    if (coord->axis != SIM_STIMULUS_AXIS_X && coord->axis != SIM_STIMULUS_AXIS_Y) {
        coord->axis = SIM_STIMULUS_AXIS_X;
    }
    if (coord->combine != SIM_STIMULUS_SEPARABLE_MULTIPLY &&
        coord->combine != SIM_STIMULUS_SEPARABLE_ADD) {
        coord->combine = SIM_STIMULUS_SEPARABLE_MULTIPLY;
    }
    if (!isfinite(coord->angle)) {
        coord->angle = 0.0;
    }
    if (!isfinite(coord->origin_x)) {
        coord->origin_x = 0.0;
    }
    if (!isfinite(coord->origin_y)) {
        coord->origin_y = 0.0;
    }
    if (!isfinite(coord->spacing_x) || fabs(coord->spacing_x) <= STIMULUS_COORD_EPS) {
        coord->spacing_x = 1.0;
    }
    if (!isfinite(coord->spacing_y) || fabs(coord->spacing_y) <= STIMULUS_COORD_EPS) {
        coord->spacing_y = coord->spacing_x;
    }
    if (!isfinite(coord->center_x)) {
        coord->center_x = 0.0;
    }
    if (!isfinite(coord->center_y)) {
        coord->center_y = 0.0;
    }
    if (!isfinite(coord->velocity_x)) {
        coord->velocity_x = 0.0;
    }
    if (!isfinite(coord->velocity_y)) {
        coord->velocity_y = 0.0;
    }
    if (!isfinite(coord->ellipse_u) || fabs(coord->ellipse_u) <= STIMULUS_COORD_EPS) {
        coord->ellipse_u = 1.0;
    }
    coord->ellipse_u = fabs(coord->ellipse_u);
    if (!isfinite(coord->ellipse_v) || fabs(coord->ellipse_v) <= STIMULUS_COORD_EPS) {
        coord->ellipse_v = coord->ellipse_u;
    }
    coord->ellipse_v = fabs(coord->ellipse_v);
    if (!isfinite(coord->spiral_arms)) {
        coord->spiral_arms = 1.0;
    }
    if (!isfinite(coord->spiral_pitch)) {
        coord->spiral_pitch = 1.0;
    }
    if (!isfinite(coord->spiral_phase)) {
        coord->spiral_phase = 0.0;
    }
    if (!isfinite(coord->spiral_angular_velocity)) {
        coord->spiral_angular_velocity = 0.0;
    }
    if (coord->mode == SIM_STIMULUS_COORD_SPIRAL &&
        fabs(coord->spiral_arms) <= STIMULUS_COORD_EPS &&
        fabs(coord->spiral_pitch) <= STIMULUS_COORD_EPS) {
        coord->spiral_arms = 1.0;
        coord->spiral_pitch = 1.0;
    }
}
 
static inline bool sim_stimulus_coord_is_time_invariant(const SimStimulusCoordConfig *coord) {
    if (coord == NULL) {
        return true;
    }
    if (fabs(coord->velocity_x) > STIMULUS_COORD_EPS ||
        fabs(coord->velocity_y) > STIMULUS_COORD_EPS) {
        return false;
    }
    if (coord->mode == SIM_STIMULUS_COORD_SPIRAL &&
        fabs(coord->spiral_angular_velocity) > STIMULUS_COORD_EPS) {
        return false;
    }
    return true;
}
 
static inline SimResult sim_stimulus_coord_xy_at_indices(const SimStimulusCoordConfig *coord,
                                                         size_t ix, size_t iy, double *out_x,
                                                         double *out_y) {
    if (coord == NULL || out_x == NULL || out_y == NULL) {
        return SIM_RESULT_INVALID_ARGUMENT;
    }
 
    *out_x = coord->origin_x + (double)ix * coord->spacing_x;
    *out_y = coord->origin_y + (double)iy * coord->spacing_y;
    return SIM_RESULT_OK;
}
 
static inline SimResult sim_stimulus_coord_patch_row(const SimStimulusCoordConfig *coord,
                                                     const SimFieldPatch *patch, size_t row_offset,
                                                     double t, SimStimulusCoordRow *out_row) {
    SimStimulusCoordRow row = {0};
 
    if (coord == NULL || !sim_field_patch_is_valid(patch) || out_row == NULL ||
        row_offset >= patch->height) {
        return SIM_RESULT_INVALID_ARGUMENT;
    }
 
    row.x0 = patch->x0;
    row.y0 = patch->y0 + row_offset;
    row.width = patch->width;
    row.x_step = coord->spacing_x;
    row.y_step = 0.0;
    row.sample_x_step = row.x_step;
    row.sample_y_step = row.y_step;
 
    if (sim_stimulus_coord_xy_at_indices(coord, row.x0, row.y0, &row.x, &row.y) != SIM_RESULT_OK) {
        return SIM_RESULT_INVALID_ARGUMENT;
    }
    sim_stimulus_coord_sample_xy(coord, row.x, row.y, t, &row.sample_x, &row.sample_y);
    *out_row = row;
    return SIM_RESULT_OK;
}
 
static inline SimResult sim_stimulus_coord_xy(const SimStimulusCoordConfig *coord,
                                              const SimField *field, size_t index, double *out_x,
                                              double *out_y) {
    size_t ix = 0U;
    size_t iy = 0U;
    SimResult rc = SIM_RESULT_OK;
 
    if (coord == NULL || field == NULL || out_x == NULL || out_y == NULL) {
        return SIM_RESULT_INVALID_ARGUMENT;
    }
 
    rc = sim_field_index_to_xy(field, index, &ix, &iy);
    if (rc != SIM_RESULT_OK) {
        return rc;
    }
 
    return sim_stimulus_coord_xy_at_indices(coord, ix, iy, out_x, out_y);
}
 
static inline double sim_stimulus_coord_u(const SimStimulusCoordConfig *coord, double x, double y,
                                          double t) {
    if (coord == NULL) {
        return x;
    }
    double sample_x = x;
    double sample_y = y;
    sim_stimulus_coord_sample_xy(coord, x, y, t, &sample_x, &sample_y);
    switch (coord->mode) {
    case SIM_STIMULUS_COORD_AXIS:
        return (coord->axis == SIM_STIMULUS_AXIS_Y) ? sample_y : sample_x;
    case SIM_STIMULUS_COORD_ANGLE: {
        double u = 0.0;
        double v = 0.0;
        sim_stimulus_coord_rotate_xy(sample_x, sample_y, coord->angle, &u, &v);
        return u;
    }
    case SIM_STIMULUS_COORD_RADIAL: {
        double dx = 0.0;
        double dy = 0.0;
        sim_stimulus_coord_centered_xy(coord, x, y, t, &dx, &dy);
        return hypot(dx, dy);
    }
    case SIM_STIMULUS_COORD_POLAR: {
        double r = 0.0;
        sim_stimulus_coord_polar(coord, x, y, t, &r, NULL);
        return r;
    }
    case SIM_STIMULUS_COORD_AZIMUTH: {
        double dx = 0.0;
        double dy = 0.0;
        sim_stimulus_coord_centered_xy(coord, x, y, t, &dx, &dy);
        return atan2(dy, dx);
    }
    case SIM_STIMULUS_COORD_ELLIPTIC: {
        double dx = 0.0;
        double dy = 0.0;
        double r = 0.0;
        sim_stimulus_coord_centered_xy(coord, x, y, t, &dx, &dy);
        sim_stimulus_coord_elliptic_polar(coord, dx, dy, &r, NULL);
        return r;
    }
    case SIM_STIMULUS_COORD_SPIRAL: {
        double dx = 0.0;
        double dy = 0.0;
        sim_stimulus_coord_centered_xy(coord, x, y, t, &dx, &dy);
        double r = hypot(dx, dy);
        double theta = atan2(dy, dx);
        return coord->spiral_pitch * r + coord->spiral_arms * theta + coord->spiral_phase +
               coord->spiral_angular_velocity * t;
    }
    case SIM_STIMULUS_COORD_SEPARABLE:
    default:
        return sample_x;
    }
}
 
#ifdef __cplusplus
}
#endif
 
#endif /* OAKFIELD_STIMULUS_COORDS_H */