/*
* stratus.c
*
* Copyright (C) 2019, Tristan Dagobert, CMLA, École Normale Supérieure Paris-Saclay.
*
* This software is a computer program.
* It is a part of a code which estimates, given a satellite time series
* the visibility of each image. Given a gray level time series, it provides a
* Boolean times series indicating the visibility (mainly the opaque clouds,
* haze and shadow).
*
*
* This software is governed by the CeCILL-C license under French law and
* abiding by the rules of distribution of free software.  You can  use,
* modify and/ or redistribute the software under the terms of the CeCILL-C
* license as circulated by CEA, CNRS and INRIA at the following URL
* "http://www.cecill.info".
*
* As a counterpart to the access to the source code and  rights to copy,
* modify and redistribute granted by the license, users are provided only
* with a limited warranty  and the software's author,  the holder of the
* economic rights,  and the successive licensors  have only  limited
* liability.
*
* In this respect, the user's attention is drawn to the risks associated
* with loading,  using,  modifying and/or developing or reproducing the
* software by the user in light of its specific status of free software,
* that may mean  that it is complicated to manipulate,  and  that  also
* therefore means  that it is reserved for developers  and  experienced
* professionals having in-depth computer knowledge. Users are therefore
* encouraged to load and test the software's suitability as regards their
* requirements in conditions enabling the security of their systems and/or
* data to be ensured and,  more generally, to use and operate it in the
* same conditions as regards security.
*
* The fact that you are presently reading this means that you have had
* knowledge of the CeCILL-C license and that you accept its terms.
*
*/

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <float.h>
#include <assert.h>

#include "iio.h"
#include "refine.h"
#include "stratus.h"
#include "utilities.h"
#include "criterion_sift.h"

/*============================================================================*/
static void
initialize_mappes_t(mappes_t *m)
{
    m->sift = NULL;
    return;
}
/*============================================================================*/
static void
create_and_link_maps_to_criterion(criterion_t * head_cn)
{
    int i, n, npix;

    criterion_t *cur;
    float **maps;

    cur = head_cn;
    npix = cur->nrow * cur->ncol;
    maps = malloc (cur->nb_ims * sizeof (float *));

    for (n = 0; n < cur->nb_ims; n++)
    {
        maps[n] = malloc (npix * sizeof (float));
        for (i = 0; i < npix; i++)    maps[n][i] = IS_GROUND;
    }

    while(cur != NULL)
    {
        cur->maps = maps;
        cur = cur->next;
    }

    return;
}
/*============================================================================*/
criterion_t *
get_criterion(criterion_t * cn, int type)
{
    criterion_t *cur;

    cur = cn;

    while(cur != NULL && cur->type != type) cur = cur->next;

    if (cur == NULL)
    {
        return NULL;
    }
    else
    {
        return cur;
    }
}
/*============================================================================*/
static void
binarize_maps(float** maps, float ** probas, int nb_ims, int nrow, int ncol,
              float epsilon)
{
    int n, j, npix;

    npix = nrow * ncol;
    for (n = 0; n < nb_ims; n++)
    {
        /* threshold */
        for (j = 0; j < npix; j++)
            maps[n][j] = (probas[n][j] < epsilon) ? IS_GROUND : IS_CLOUD;

    }
    return;
}
/*============================================================================*/
/* compute the probabilities of the image series within a patches and
 * multi-shifted approach
 *============================================================================*/
mappes_t
compute_cloud_maps (criterion_t * head_cn, options_t opt)
{

    int nb_ims, nrow, ncol;
    criterion_t *cn = head_cn;

    nb_ims =cn->nb_ims;
    nrow = cn->nrow;
    ncol = cn->ncol;

    mappes_t mappes;
    initialize_mappes_t(&mappes);

    /* Boolean maps indicating, per pixel, if there a cloud or not. The boolean
     * series is linked so to be common to all criterion structures */
    create_and_link_maps_to_criterion(head_cn);

    /* compute the distances between the SIFT descriptors */
    int slig, scol;
    get_sift_image_size(&slig, &scol, nrow, ncol);

    float ** ims_corr = NULL;
    mappes.sift = alloc_series_f(nb_ims, nrow * ncol);

    while (cn != NULL)
    {
        if (cn->type != SIFT)
        {
            cn = cn->next;
            continue;
        }

        ims_corr =
            compute_distance_sift_series (cn->ims, nb_ims, nrow, ncol);
        /* binarize */
        binarize_maps(cn->maps, ims_corr, nb_ims, nrow, ncol, opt.alpha);

        free_series_f (ims_corr, nb_ims, NULL);

        ims_corr = refine_boolean_series(cn->maps, nb_ims, nrow, ncol, opt.theta);

        memcpy_series(mappes.sift, ims_corr, nb_ims, nrow * ncol);
        memcpy_series(cn->maps, ims_corr, nb_ims, nrow * ncol);
        free_series_f (ims_corr, nb_ims, NULL);

        break;
    }

    return mappes;
}