El software del detector: clase y programas

Accesibilidad: Alt+P para escuchar o pausar, Alt+S para detener.

Proyecto original: Cámara de rayos cósmicos con móvil
Autores: Francisco Albiol y Alberto Corbi — IFIC CSIC-UV
Premio: CPAN (Centro Nacional de Física de Partículas, Astropartículas y Nuclear)
Licencia: código publicado con fines educativos y de divulgación científica

El detector se implementa en C++ sobre dos librerías multiplataforma: OpenCV para la captura desde cámara y CImg para el análisis de imagen píxel a píxel. Consulta Frameworks y herramientas para instalarlas en tu sistema.

La clase `ccdRadiationDetector`

La clase encapsula toda la lógica del detector: gestión de calibraciones, adquisición de frames, identificación de hits y construcción del histograma de energía depositada.

class ccdRadiationDetector {
public:
    typedef std::map<int,int> mapcounter_t;
private:
    // Bases de datos de calibración
    cimg_library::CImg<float> *_histogram;
    cimg_library::CImg<float> *_frame;
    cimg_library::CImg<float> *_calibration;
    cimg_library::CImg<float> *_noise;
    cimg_library::CImg<float> *_badchannels;
    // Parámetros estadísticos
    size_t all_events;
    size_t hits;
    // Nombres de ficheros para persistencia de calibración
    std::string _histogram_filename;
    std::string _bighit_filename;
    std::string _maskhit_filename;
    std::string _calibration_filename;

    mapcounter_t _badchannel_counter;
    double _maxval;

    static int mask_image(const cimg_library::CImg<float> &noise,
                          const cimg_library::CImg<float> &badchannels,
                          cimg_library::CImg<float> &frame,
                          float noise_sigmas,
                          mapcounter_t &badchannel_counter);
public:
    ccdRadiationDetector(int device, const std::string &outputpath="");
    ~ccdRadiationDetector();
    bool has_calibration() const;
    void calibrate(int nframes=2048);
    void adquire(float threshold, int nframes=200);
    float modo_proporcional(float threshold, int nframes=200);
    void show_histogram(int scale=1) const;
};

Calibración del sensor

La calibración acumula nframes imágenes oscuras para calcular píxel a píxel la media (pedestales) y la desviación estándar (ruido). Los canales con relación señal/ruido anormal se marcan como badchannels y se excluyen del análisis.

void ccdRadiationDetector::calibrate(int nframes)
{
    int loopcount = nframes < 20 ? 20 : nframes;
    cimg_library::CImg<> animage, noise2, value, gry;

    _calibration->assign();
    value.load_video("Default");     // Primera imagen
    noise2 = value.get_sqr();         // Acumulación al cuadrado

    for (int i = 1; i < loopcount; i++) {
        animage.load_video("Default");
        value  += animage;
        noise2 += animage.get_sqr();
    }

    (*_calibration) = (value) * (1. / float(loopcount));  // Media (pedestales)
    {
        _badchannels->assign(*_calibration);
        _badchannels->fill(0);
        // Desviación estándar píxel a píxel
        (*_noise) = (noise2 - (value.mul(*_calibration))) / (loopcount - 1);
        _noise->sqrt();

        cimg_library::CImg<> stats      = _calibration->get_stats();
        cimg_library::CImg<> noisestats = _noise->get_stats();

        const float variance_value = sqrt(stats(4));
        float noise_level = (stats(1) - stats(0)) / variance_value;

        // Marcar píxeles ruidosos como badchannels
        cimg_foroff(*_badchannels, off) {
            const float val   = (*_calibration)(off);
            const float noise = (*_noise)(off);
            if (noise == 0 || (val / noise) > noise_level)
                (*_badchannels)(off) = 1;
        }
    }

    _badchannels->dilate(6, 6);   // Extender badchannels 6 píxeles

    // Propagar badchannels a los tres canales RGB
    cimg_forXY(*_badchannels, x, y) {
        if ((*_badchannels)(x,y,0) || (*_badchannels)(x,y,1) || (*_badchannels)(x,y,2)) {
            for (int c = 0; c < 3; c++) {
                (*_calibration)(x,y,c)  = 0;
                (*_noise)(x,y,c)        = 0;
                (*_badchannels)(x,y,c)  = 1;
            }
        }
    }

    // Mostrar resultados de calibración
    _calibration->display("Máscara de calibración (pedestales)");
    _noise->display("Ruido de los canales");
    _badchannels->display("Bad channels");

    // Guardar calibración para reutilizarla
    cimg_library::CImgList<float> calibration_list;
    calibration_list.push_back(*_calibration);
    calibration_list.push_back(*_noise);
    calibration_list.push_back(*_badchannels);
    calibration_list.save(_calibration_filename.c_str());
}

Programa de captura

Todos los programas se ejecutan en un directorio. Requieren OpenCV y CImg instalados.

//  main_capture.cpp
//  Creado por Kiko Albiol — IFIC CSIC, 13/09/14
//  Captura continua y análisis de radiación

#include <iostream>
#include "opencvImageCapture.h"

int main(int argc, const char * argv[])
{
    int device = 0;
    if (argc > 1)
        device = atoi(argv[1]);

    std::cout << "Abriendo captura de video " << device << std::endl;
    ccdRadiationDetector *rad_detector = new ccdRadiationDetector(device);

    if (rad_detector->has_calibration() == false)
        rad_detector->calibrate(800);  // 800 imágenes para calcular ruido y pedestales

    while (1) {
        rad_detector->adquire(1000);   // Adquiere 1000 imágenes, procesa y
                                        // si encuentra señal la guarda en el histograma
    }
}

Programa de visualización periódica del histograma

//  main_display.cpp
//  Creado por Kiko Albiol — IFIC CSIC, 13/09/14
//  Muestra el histograma de energía acumulada

#include <iostream>
#include "opencvImageCapture.h"

int main(int argc, const char * argv[])
{
    ccdRadiationDetector *rad_detector = new ccdRadiationDetector(0);

    while (1) {
        // Comprime el histograma en 16 partes para visualización
        rad_detector->show_histogram(16);
    }
}

Visor de imagen simple

Útil para inspeccionar visualmente los ficheros .cimg generados por la calibración.

//  main_viewer.cpp
//  Creado por Kiko Albiol — IFIC CSIC, 13/09/14
//  Visualizador de imágenes CImg

#include <iostream>
#include "CImg.h"

int main(int argc, const char * argv[])
{
    if (argc < 2) {
        std::cout << "[" << argv[0] << "]  imagen.cimg" << std::endl;
        return 1;
    }
    const std::string filename(argv[1]);
    cimg_library::CImg<> animage(filename.c_str());
    animage.display("Image Displayed...");
}