CLSMImage.h

Go to the documentation of this file.

#ifndef TTTRLIB_CLSMIMAGE_H
#define TTTRLIB_CLSMIMAGE_H
 
#include <iostream> /* cout, clog */
#include <vector>
#include <utility>
 
#include <cstring>
#include <complex>
#include <cmath>
 
#include "pocketfft/pocketfft_hdronly.h"
 
#include "TTTR.h" /* TTTR */
 
#include "CLSMFrame.h"
#include "CLSMLine.h"
#include "CLSMPixel.h"
#include "DecayPhasor.h"
#include "Correlator.h"
 
 
typedef enum{
    CLSM_DEFAULT,         
    CLSM_SP5,             
    CLSM_SP8              
} ReadingRoutine;
 
 
template<typename TS, typename TE>
static std::pair<int, int> find_clsm_start_stop(
        int &i_event,
        int marker_start, int marker_stop, int marker_event,
        unsigned long long* macro_time_arr,
        TS start_stop_arr,
        TE event_type_arr,
        size_t n,
        long duration = -1
){
    int start = -1;
    int stop = -1;
    for(; i_event < n; i_event++)
    {
        if(event_type_arr[i_event] != marker_event) continue;
        if(start_stop_arr[i_event] == marker_start)
            // start found -> search stop and break
        {
            start = i_event;
            if(marker_stop > 0) // search for stop idx using stop event
            {
                for(;i_event < n; i_event++){
                    if(event_type_arr[i_event] != marker_event) continue;
                    if(start_stop_arr[i_event] == marker_stop){
                        stop = i_event;
                        break;
                    }
                }
            }
            else if(duration > 0) // search for stop idx using duration
            {
                unsigned long stop_time = macro_time_arr[i_event] + duration;
                for(;i_event < n; i_event++){
                    if(macro_time_arr[i_event] >= stop_time){
                        stop = i_event;
                        break;
                    }
                }
            }
            else // do not search for stop idx
            {
                stop = 0;
            }
            break;
        }
    }
    return {start, stop};
}
 
 
 
class CLSMSettings{
 
    friend class CLSMImage;
 
protected:
 
    bool skip_before_first_frame_marker = false;
    bool skip_after_last_frame_marker = false;
 
    int reading_routine = CLSM_DEFAULT;
 
    int marker_line_start = 0;
 
    int marker_line_stop = 0;
 
    std::vector<int> marker_frame_start = {};
 
    int marker_event_type = 0;
 
    int n_pixel_per_line = 0;
    int n_lines = 0;
 
public:
 
    explicit CLSMSettings(
        bool skip_before_first_frame_marker = false,
        bool skip_after_last_frame_marker = false,
        int reading_routine = CLSM_DEFAULT,
        int marker_line_start = 3,
        int marker_line_stop = 2,
        std::vector<int> marker_frame_start = std::vector<int>({1}),
        int marker_event_type = 1,
        int n_pixel_per_line = 1,
        int n_lines = -1
        // long long macro_time_shift = 0
    ){
        this->skip_before_first_frame_marker = skip_before_first_frame_marker;
        this->skip_after_last_frame_marker = skip_after_last_frame_marker;
        this->reading_routine = reading_routine;
        this->n_pixel_per_line = n_pixel_per_line;
        this->n_lines = n_lines;
        this->marker_event_type = marker_event_type;
        this->marker_line_stop = marker_line_stop;
        this->marker_line_start = marker_line_start;
        this->marker_frame_start = marker_frame_start;
//        this->macro_time_shift = macro_time_shift;
    }
 
};
 
 
class CLSMImage {
 
    friend class Correlator;
    friend class CLSMFrame;
    friend class CLSMLine;
    friend class CLSMPixel;
 
private:
 
    CLSMSettings settings;
 
    bool _is_filled_ = false;
 
    std::vector<CLSMFrame *> frames;
 
    void remove_incomplete_frames();
 
    void create_pixels_in_lines();
 
protected:
 
    size_t n_frames = 0;
 
    size_t n_lines = 0;
 
    size_t n_pixel = 0;
 
    std::shared_ptr<TTTR> tttr = nullptr;
 
    void create_frames(bool clear_first = true);
 
    void create_lines();
 
    void determine_number_of_lines();
 
public:
 
    std::shared_ptr<TTTR> get_tttr(){
        return tttr;
    }
 
    void set_tttr(std::shared_ptr<TTTR> v){
        tttr = v;
    }    
 
    const CLSMSettings* get_settings(){
        return &settings;
    }
 
    size_t size(){
        return frames.size();
    }
 
    std::vector<int> marker_frame;
 
    int marker_line_start;
 
    int marker_line_stop;
 
    int marker_event;
 
    std::string reading_routine = "default";
 
    bool skip_before_first_frame_marker = false;
    bool skip_after_last_frame_marker = false;
 
    void fill(
            TTTR *tttr_data = nullptr,
            std::vector<int> channels = std::vector<int>(),
            bool clear = true,
            const std::vector<std::pair<int,int>> &micro_time_ranges = std::vector<std::pair<int,int>>()
    );
 
    void fill_pixels(
            TTTR *tttr_data,
            std::vector<int> channels,
            bool clear_pixel = true,
            std::vector<std::pair<int,int>> micro_time_ranges = std::vector<std::pair<int,int>>()
    ){
        std::clog << "WARNING: 'fill_pixels' deprecated.  Use 'fill'." << std::endl;
        fill(tttr_data, channels, clear_pixel, micro_time_ranges);
    }
 
    void clear();
 
    void clear_pixels(){
        std::clog << "WARNING: 'clear_pixels' deprecated.  Use 'clear'." << std::endl;
        clear();
    }
 
    void strip(const std::vector<int> &tttr_indices, int offset = 0);
 
    void get_fcs_image(
        float **output, int *dim1, int *dim2, int *dim3, int *dim4,
        std::shared_ptr<TTTR> tttr,
        CLSMImage *clsm_other = nullptr,
        std::string correlation_method = "default",
        int n_bins = 50,
        int n_casc = 1,
        bool stack_frames = false,
        bool normalized_correlation = false,
        int min_photons = 2
    );
 
    std::vector<CLSMFrame *> get_frames() {
        return frames;
    }
 
    void get_intensity(unsigned short **output, int *dim1, int *dim2, int *dim3);
 
    void get_fluorescence_decay(
        TTTR *tttr_data,
        unsigned char **output, int *dim1, int *dim2, int *dim3, int *dim4,
        int micro_time_coarsening = 1,
        bool stack_frames = false
    );
 
 
    void get_decay_of_pixels(
        TTTR *tttr_data,
        uint8_t* mask, int dmask1, int dmask2, int dmask3,
        unsigned int **output, int *dim1, int *dim2,
        int tac_coarsening,
        bool stack_frames
    );
 
 
    void get_mean_micro_time(
        TTTR *tttr_data,
        double **output, int *dim1, int *dim2, int *dim3,
        double microtime_resolution = -1.0,
        int minimum_number_of_photons = 2,
        bool stack_frames = false
    );
 
 
    void get_phasor(
        float **output, int *dim1, int *dim2, int *dim3, int *dim4,
        TTTR *tttr_data,
        TTTR *tttr_irf = nullptr,
        double frequency = -1,
        int minimum_number_of_photons = 2,
        bool stack_frames = false
    );
 
 
    void get_mean_lifetime(
            TTTR *tttr_data,
            double **output, int *dim1, int *dim2, int *dim3,
            int minimum_number_of_photons = 3,
            TTTR *tttr_irf = nullptr, double m0_irf = 1.0, double m1_irf = 1.0,
            bool stack_frames = false,
            std::vector<double> background = std::vector<double>(),
            double m0_bg = 0.0, double m1_bg = 0.0,
            double background_fraction = -1.0
    );
 
    inline int to1D(int frame, int line, int pixel) {
        return (frame * n_lines * n_pixel) + (line * n_pixel) + pixel;
    }
 
 
    inline std::vector<int> to3D(int idx) {
        int frame = idx / (n_lines * n_pixel);
        idx -= (frame * n_lines * n_pixel);
        int line = idx / n_pixel;
        int pixel = idx % n_pixel;
        return std::vector<int>{frame, line, pixel};
    }
 
 
    CLSMPixel* getPixel(unsigned int idx) {
        int frame, line, pixel;
 
        frame = idx / (n_lines * n_pixel);
        idx -= (frame * n_lines * n_pixel);
        line = idx / n_pixel;
        pixel = idx % n_pixel;
 
        CLSMFrame* s_frame = frames[frame];
        CLSMLine*  s_line  = s_frame->lines[line];
        CLSMPixel* s_pixel = &(s_line->pixels[pixel]);
        return s_pixel;
    }
 
 
    int get_n_frames() const {
        return n_frames;
    }
 
 
    int get_n_lines() const {
        return n_lines;
    }
 
 
    int get_n_pixel() const {
        return n_pixel;
    }
 
    void copy(const CLSMImage &p2, bool fill = false);
 
    void append(CLSMFrame *frame);
 
    void transform(unsigned int* input, int n_input);
 
    void rebin(int bin_line, int bin_pixel);
 
    void distribute(
        unsigned int pixel_id,
        CLSMImage* target,
        std::vector<int>& target_pixel_ids,
        std::vector<int>& target_probabilities
    );
 
 
    void crop(
        int frame_start, int frame_stop,
        int line_start, int line_stop,
        int pixel_start, int pixel_stop
    );
 
    void stack_frames() {
        CLSMFrame* f0 = frames[0];
        for (unsigned int i = 1; i < n_frames; i++) {
            *f0 += *frames[i];
        }
        frames.resize(1);
        n_frames = 1;
    }
 
 
    CLSMImage(const CLSMImage &p2, bool fill = false);
 
    explicit CLSMImage(
            std::shared_ptr<TTTR> tttr_data = nullptr,
            CLSMSettings settings = CLSMSettings(),
            CLSMImage *source = nullptr,
            bool fill = true,
            std::vector<int> channels = std::vector<int>(),
            std::vector<std::pair<int, int>> micro_time_ranges =
                    std::vector<std::pair<int, int>>()
    );
 
    virtual ~CLSMImage() {
        for (auto frame : frames) {
            delete frame;
        }
    }
 
 
    CLSMFrame* operator[](unsigned int i_frame) {
        return frames[i_frame];
    }
 
 
    static void compute_ics(
            double **output, int *dim1, int *dim2, int *dim3,
            std::shared_ptr<TTTR> tttr_data = nullptr,
            CLSMImage* clsm = nullptr,
            double *images = nullptr, int input_frames = -1, int input_lines = -1, int input_pixel = 1,
            std::vector<int> x_range = std::vector<int>({0, -1}),
            std::vector<int> y_range = std::vector<int>({0, -1}),
            std::vector<std::pair<int, int>> frames_index_pairs = std::vector<std::pair<int, int>>(),
            std::string subtract_average = "",
            uint8_t *mask = nullptr, int dmask1 = -1, int dmask2 = -1, int dmask3 = -1
    );
 
    static void get_roi(
             double** output, int* dim1, int* dim2, int* dim3,
             CLSMImage* clsm = nullptr,
             std::vector<int> x_range=std::vector<int>({0,-1}),
             std::vector<int> y_range=std::vector<int>({0,-1}),
             std::string subtract_average = "",
             double background = 0.0,
             bool clip=false, double clip_max=1e6, double clip_min=-1e6,
             double *images = nullptr, int n_frames=-1, int n_lines=-1, int n_pixels=1,
             uint8_t *mask = nullptr, int dmask1 = -1, int dmask2 = -1, int dmask3 = -1,
             std::vector<int> selected_frames = std::vector<int>()
    );
 
 
    static std::vector<int> get_frame_edges(
            TTTR* tttr = nullptr,
            int start_event = 0,
            int stop_event = -1,
            std::vector<int> marker_frame_start = std::vector<int>({4, 6}),
            int marker_event_type = 15,
            int reading_routine = CLSM_SP8,
            bool skip_before_first_frame_marker = false,
            bool skip_after_last_frame_marker = false
    );
 
    static std::vector<int> get_line_edges(
        TTTR* tttr,
        int start_event, int stop_event,
        int marker_line_start = 1, int marker_line_stop = 2,
        int marker_event_type = 15,
        int reading_routine = CLSM_SP8
    );
 
 
    static std::vector<int> get_line_edges_by_duration(
        TTTR* tttr,
        int start_event, int stop_event,
        int marker_line_start = 1,
        int line_duration = 2,
        int marker_event_type = 15,
        int reading_routine = CLSM_SP8
    );
 
 
    double get_line_duration(int frame = 0, int line = 0){
        double re = -1.0;
        if(tttr != nullptr){
            auto header = tttr->get_header();
            
            auto f = frames[frame];
            auto l = f->lines[line];
 
            int start = l->get_start();
            int stop = l->get_stop();
 
            unsigned long long t_stop = tttr->macro_times[stop];
            unsigned long long t_start = tttr->macro_times[start];
            unsigned long long dt = t_stop - t_start;
            double res = header->get_macro_time_resolution() * 1000.0;
            re = dt * res;
        }
        return re;
    }
 
    double get_pixel_duration(int frame = 0, int line = 0){
        return get_line_duration(frame, line) / settings.n_pixel_per_line;
    }
 
};
 
 
#endif //TTTRLIB_CLSMIMAGE_H