cvaux/src/cvhmm1d.cpp

/*M///////////////////////////////////////////////////////////////////////////////////////
//
//  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
//
//  By downloading, copying, installing or using the software you agree to this license.
//  If you do not agree to this license, do not download, install,
//  copy or use the software.
//
//
//                        Intel License Agreement
//                For Open Source Computer Vision Library
//
// Copyright (C) 2000, Intel Corporation, all rights reserved.
// Third party copyrights are property of their respective owners.
//
// Redistribution and use in source and binary forms, with or without modification,
// are permitted provided that the following conditions are met:
//
//   * Redistribution's of source code must retain the above copyright notice,
//     this list of conditions and the following disclaimer.
//
//   * Redistribution's in binary form must reproduce the above copyright notice,
//     this list of conditions and the following disclaimer in the documentation
//     and/or other materials provided with the distribution.
//
//   * The name of Intel Corporation may not be used to endorse or promote products
//     derived from this software without specific prior written permission.
//
// This software is provided by the copyright holders and contributors "as is" and
// any express or implied warranties, including, but not limited to, the implied
// warranties of merchantability and fitness for a particular purpose are disclaimed.
// In no event shall the Intel Corporation or contributors be liable for any direct,
// indirect, incidental, special, exemplary, or consequential damages
// (including, but not limited to, procurement of substitute goods or services;
// loss of use, data, or profits; or business interruption) however caused
// and on any theory of liability, whether in contract, strict liability,
// or tort (including negligence or otherwise) arising in any way out of
// the use of this software, even if advised of the possibility of such damage.
//
//M*/


#include "_cvaux.h"

#if 0

#define LN2PI 1.837877f
#define BIG_FLT 1.e+10f


#define _CV_ERGODIC 1
#define _CV_CAUSAL 2

#define _CV_LAST_STATE 1
#define _CV_BEST_STATE 2

//*F///////////////////////////////////////////////////////////////////////////////////////
//    Name: icvForward1DHMM
//    Purpose: The function performs baum-welsh algorithm
//    Context:
//    Parameters: obs_info - addres of pointer to CvImgObsInfo structure
//                num_hor_obs - number of horizontal observation vectors
//                num_ver_obs - number of horizontal observation vectors
//                obs_size - length of observation vector
//
//    Returns: error status
//
//    Notes:
//F*/
#if 0
CvStatus icvForward1DHMM( int num_states, int num_obs, CvMatr64d A,
                          CvMatr64d B,
                          double* scales)
{
    // assume that observation and transition
    // probabilities already computed
    int m_HMMType  = _CV_CAUSAL;
    double* m_pi = icvAlloc( num_states* sizeof( double) );

    /* alpha is matrix
       rows throuhg states
       columns through time
    */
    double* alpha = icvAlloc( num_states*num_obs * sizeof( double ) );

    /* All calculations will be in non-logarithmic domain */

    /* Initialization */
    /* set initial state probabilities */
    m_pi[0] = 1;
    for (i = 1; i < num_states; i++)
    {
        m_pi[i] = 0.0;
    }

    for  (i = 0; i < num_states; i++)
    {
        alpha[i] = m_pi[i] * m_b[ i];
    }

    /******************************************************************/
    /*   Induction                                                    */

    if ( m_HMMType == _CV_ERGODIC )
    {
        int t;
        for (t = 1 ; t < num_obs; t++)
        {
            for (j = 0; j < num_states; j++)
            {
               double sum = 0.0;
               int i;

                for (i = 0; i < num_states; i++)
                {
                     sum += alpha[(t - 1) * num_states + i] * A[i * num_states + j];
                }

                alpha[(t - 1) * num_states + j] = sum * B[t * num_states + j];

                /* add computed alpha to scale factor */
                sum_alpha += alpha[(t - 1) * num_states + j];
            }

            double scale = 1/sum_alpha;

            /* scale alpha */
            for (j = 0; j < num_states; j++)
            {
                alpha[(t - 1) * num_states + j] *= scale;
            }

            scales[t] = scale;

        }
    }

#endif


//*F///////////////////////////////////////////////////////////////////////////////////////
//    Name: icvCreateObsInfo
//    Purpose: The function allocates memory for CvImgObsInfo structure
//             and its inner stuff
//    Context:
//    Parameters: obs_info - addres of pointer to CvImgObsInfo structure
//                num_hor_obs - number of horizontal observation vectors
//                num_ver_obs - number of horizontal observation vectors
//                obs_size - length of observation vector
//
//    Returns: error status
//
//    Notes:
//F*/
/*CvStatus icvCreateObsInfo( CvImgObsInfo** obs_info,
                              CvSize num_obs, int obs_size )
{
    int total = num_obs.height * num_obs.width;

    CvImgObsInfo* obs = (CvImgObsInfo*)icvAlloc( sizeof( CvImgObsInfo) );

    obs->obs_x = num_obs.width;
    obs->obs_y = num_obs.height;

    obs->obs = (float*)icvAlloc( total * obs_size * sizeof(float) );

    obs->state = (int*)icvAlloc( 2 * total * sizeof(int) );
    obs->mix = (int*)icvAlloc( total * sizeof(int) );

    obs->obs_size = obs_size;

    obs_info[0] = obs;

    return CV_NO_ERR;
}*/

/*CvStatus icvReleaseObsInfo( CvImgObsInfo** p_obs_info )
{
    CvImgObsInfo* obs_info = p_obs_info[0];

    icvFree( &(obs_info->obs) );
    icvFree( &(obs_info->mix) );
    icvFree( &(obs_info->state) );
    icvFree( &(obs_info) );

    p_obs_info[0] = NULL;

    return CV_NO_ERR;
} */


//*F///////////////////////////////////////////////////////////////////////////////////////
//    Name: icvCreate1DHMM
//    Purpose: The function allocates memory for 1-dimensional HMM
//             and its inner stuff
//    Context:
//    Parameters: hmm - addres of pointer to CvEHMM structure
//                state_number - number of states in HMM
//                num_mix - number of gaussian mixtures in HMM states
//                          size of array is defined by previous parameter
//                obs_size - length of observation vectors
//
//    Returns: error status
//    Notes:
//F*/
CvStatus icvCreate1DHMM( CvEHMM** this_hmm,
                         int state_number, int* num_mix, int obs_size )
{
    int i;
    int real_states = state_number;

    CvEHMMState* all_states;
    CvEHMM* hmm;
    int total_mix = 0;
    float* pointers;

    /* allocate memory for hmm */
    hmm = (CvEHMM*)icvAlloc( sizeof(CvEHMM) );

    /* set number of superstates */
    hmm->num_states = state_number;
    hmm->level = 0;

    /* allocate memory for all states */
    all_states = (CvEHMMState *)icvAlloc( real_states * sizeof( CvEHMMState ) );

    /* assign number of mixtures */
    for( i = 0; i < real_states; i++ )
    {
        all_states[i].num_mix = num_mix[i];
    }

    /* compute size of inner of all real states */
    for( i = 0; i < real_states; i++ )
    {
        total_mix += num_mix[i];
    }
    /* allocate memory for states stuff */
    pointers = (float*)icvAlloc( total_mix * (2/*for mu invvar */ * obs_size +
                                 2/*for weight and log_var_val*/ ) * sizeof( float) );

    /* organize memory */
    for( i = 0; i < real_states; i++ )
    {
        all_states[i].mu      = pointers; pointers += num_mix[i] * obs_size;
        all_states[i].inv_var = pointers; pointers += num_mix[i] * obs_size;

        all_states[i].log_var_val = pointers; pointers += num_mix[i];
        all_states[i].weight      = pointers; pointers += num_mix[i];
    }
    hmm->u.state = all_states;

    hmm->transP = icvCreateMatrix_32f( hmm->num_states, hmm->num_states );
    hmm->obsProb = NULL;

    /* if all ok - return pointer */
    *this_hmm = hmm;
    return CV_NO_ERR;
}

CvStatus icvRelease1DHMM( CvEHMM** phmm )
{
    CvEHMM* hmm = phmm[0];
    icvDeleteMatrix( hmm->transP );

    if (hmm->obsProb != NULL)
    {
        int* tmp = ((int*)(hmm->obsProb)) - 3;
        icvFree( &(tmp)  );
    }

    icvFree( &(hmm->u.state->mu) );
    icvFree( &(hmm->u.state) );

    phmm[0] = NULL;

    return CV_NO_ERR;
}

/*can be used in CHMM & DHMM */
CvStatus icvUniform1DSegm( Cv1DObsInfo* obs_info, CvEHMM* hmm )
{
    /* implementation is very bad */
    int  i;
    CvEHMMState* first_state;

    /* check arguments */
    if ( !obs_info || !hmm ) return CV_NULLPTR_ERR;

    first_state = hmm->u.state;

    for (i = 0; i < obs_info->obs_x; i++)
    {
        //bad line (division )
        int state = (i * hmm->num_states)/obs_info->obs_x;
        obs_info->state[i] = state;
    }
    return CV_NO_ERR;
}


/*F///////////////////////////////////////////////////////////////////////////////////////
//    Name: InitMixSegm
//    Purpose: The function implements the mixture segmentation of the states of the embedded HMM
//    Context: used with the Viterbi training of the embedded HMM
//             Function uses K-Means algorithm for clustering
//
//    Parameters:  obs_info_array - array of pointers to image observations
//                 num_img - length of above array
//                 hmm - pointer to HMM structure
//
//    Returns: error status
//
//    Notes:
//F*/
CvStatus icvInit1DMixSegm(Cv1DObsInfo** obs_info_array, int num_img, CvEHMM* hmm)
{
    int  k, i, j;
    int* num_samples; /* number of observations in every state */
    int* counter;     /* array of counters for every state */

    int**  a_class;   /* for every state - characteristic array */

    CvVect32f** samples; /* for every state - pointer to observation vectors */
    int***  samples_mix;   /* for every state - array of pointers to vectors mixtures */

    CvTermCriteria criteria = cvTermCriteria( CV_TERMCRIT_EPS|CV_TERMCRIT_ITER,
                                              1000,    /* iter */
                                              0.01f ); /* eps  */

    int total = hmm->num_states;
    CvEHMMState* first_state = hmm->u.state;

    /* for every state integer is allocated - number of vectors in state */
    num_samples = (int*)icvAlloc( total * sizeof(int) );

    /* integer counter is allocated for every state */
    counter = (int*)icvAlloc( total * sizeof(int) );

    samples = (CvVect32f**)icvAlloc( total * sizeof(CvVect32f*) );
    samples_mix = (int***)icvAlloc( total * sizeof(int**) );

    /* clear */
    memset( num_samples, 0 , total*sizeof(int) );
    memset( counter, 0 , total*sizeof(int) );


    /* for every state the number of vectors which belong to it is computed (smth. like histogram) */
    for (k = 0; k < num_img; k++)
    {
        CvImgObsInfo* obs = obs_info_array[k];

        for (i = 0; i < obs->obs_x; i++)
        {
            int state = obs->state[ i ];
            num_samples[state] += 1;
        }
    }

    /* for every state int* is allocated */
    a_class = (int**)icvAlloc( total*sizeof(int*) );

    for (i = 0; i < total; i++)
    {
        a_class[i] = (int*)icvAlloc( num_samples[i] * sizeof(int) );
        samples[i] = (CvVect32f*)icvAlloc( num_samples[i] * sizeof(CvVect32f) );
        samples_mix[i] = (int**)icvAlloc( num_samples[i] * sizeof(int*) );
    }

    /* for every state vectors which belong to state are gathered */
    for (k = 0; k < num_img; k++)
    {
        CvImgObsInfo* obs = obs_info_array[k];
        int num_obs = obs->obs_x;
        float* vector = obs->obs;

        for (i = 0; i < num_obs; i++, vector+=obs->obs_size )
        {
            int state = obs->state[i];

            samples[state][counter[state]] = vector;
            samples_mix[state][counter[state]] = &(obs->mix[i]);
            counter[state]++;
        }
    }

    /* clear counters */
    memset( counter, 0, total*sizeof(int) );

    /* do the actual clustering using the K Means algorithm */
    for (i = 0; i < total; i++)
    {
        if ( first_state[i].num_mix == 1)
        {
            for (k = 0; k < num_samples[i]; k++)
            {
                /* all vectors belong to one mixture */
                a_class[i][k] = 0;
            }
        }
        else if( num_samples[i] )
        {
            /* clusterize vectors  */
            icvKMeans( first_state[i].num_mix, samples[i], num_samples[i],
                obs_info_array[0]->obs_size, criteria, a_class[i] );
        }
    }

    /* for every vector number of mixture is assigned */
    for( i = 0; i < total; i++ )
    {
        for (j = 0; j < num_samples[i]; j++)
        {
            samples_mix[i][j][0] = a_class[i][j];
        }
    }

   for (i = 0; i < total; i++)
    {
        icvFree( &(a_class[i]) );
        icvFree( &(samples[i]) );
        icvFree( &(samples_mix[i]) );
    }

    icvFree( &a_class );
    icvFree( &samples );
    icvFree( &samples_mix );
    icvFree( &counter );
    icvFree( &num_samples );


    return CV_NO_ERR;
}

/*F///////////////////////////////////////////////////////////////////////////////////////
//    Name: ComputeUniModeGauss
//    Purpose: The function computes the Gaussian pdf for a sample vector
//    Context:
//    Parameters:  obsVeq - pointer to the sample vector
//                 mu - pointer to the mean vector of the Gaussian pdf
//                 var - pointer to the variance vector of the Gaussian pdf
//                 VecSize - the size of sample vector
//
//    Returns: the pdf of the sample vector given the specified Gaussian
//
//    Notes:
//F*/
/*float icvComputeUniModeGauss(CvVect32f vect, CvVect32f mu,
                              CvVect32f inv_var, float log_var_val, int vect_size)
{
    int n;
    double tmp;
    double prob;

    prob = -log_var_val;

    for (n = 0; n < vect_size; n++)
    {
        tmp = (vect[n] - mu[n]) * inv_var[n];
        prob = prob - tmp * tmp;
   }
   //prob *= 0.5f;

   return (float)prob;
}*/

/*F///////////////////////////////////////////////////////////////////////////////////////
//    Name: ComputeGaussMixture
//    Purpose: The function computes the mixture Gaussian pdf of a sample vector.
//    Context:
//    Parameters:  obsVeq - pointer to the sample vector
//                 mu  - two-dimensional pointer to the mean vector of the Gaussian pdf;
//                       the first dimension is indexed over the number of mixtures and
//                       the second dimension is indexed along the size of the mean vector
//                 var - two-dimensional pointer to the variance vector of the Gaussian pdf;
//                       the first dimension is indexed over the number of mixtures and
//                       the second dimension is indexed along the size of the variance vector
//                 VecSize - the size of sample vector
//                 weight - pointer to the wights of the Gaussian mixture
//                 NumMix - the number of Gaussian mixtures
//
//    Returns: the pdf of the sample vector given the specified Gaussian mixture.
//
//    Notes:
//F*/
/* Calculate probability of observation at state in logarithmic scale*/
/*float icvComputeGaussMixture( CvVect32f vect, float* mu,
                                float* inv_var, float* log_var_val,
                                int vect_size, float* weight, int num_mix )
{
    double prob, l_prob;

    prob = 0.0f;

    if (num_mix == 1)
    {
        return icvComputeUniModeGauss( vect, mu, inv_var, log_var_val[0], vect_size);
    }
    else
    {
        int m;
        for (m = 0; m < num_mix; m++)
        {
            if ( weight[m] > 0.0)
            {
                l_prob = icvComputeUniModeGauss(vect, mu + m*vect_size,
                                                        inv_var + m * vect_size,
                                                        log_var_val[m],
                                                        vect_size);

                prob = prob + weight[m]*exp((double)l_prob);
            }
        }
        prob = log(prob);
    }
    return (float)prob;
}
*/

/*F///////////////////////////////////////////////////////////////////////////////////////
//    Name: EstimateObsProb
//    Purpose: The function computes the probability of every observation in every state
//    Context:
//    Parameters:  obs_info - observations
//                 hmm      - hmm
//    Returns: error status
//
//    Notes:
//F*/
CvStatus icvEstimate1DObsProb(CvImgObsInfo* obs_info, CvEHMM* hmm )
{
    int j;
    int total_states = 0;

    /* check if matrix exist and check current size
       if not sufficient - realloc */
    int status = 0; /* 1 - not allocated, 2 - allocated but small size,
                       3 - size is enough, but distribution is bad, 0 - all ok */

    /*for( j = 0; j < hmm->num_states; j++ )
    {
       total_states += hmm->u.ehmm[j].num_states;
    }*/
    total_states = hmm->num_states;

    if ( hmm->obsProb == NULL )
    {
        /* allocare memory */
        int need_size = ( obs_info->obs_x /* * obs_info->obs_y*/ * total_states * sizeof(float) /* +
                          obs_info->obs_y * hmm->num_states * sizeof( CvMatr32f) */);

        int* buffer = (int*)icvAlloc( need_size + 3 * sizeof(int) );
        buffer[0] = need_size;
        buffer[1] = obs_info->obs_y;
        buffer[2] = obs_info->obs_x;
        hmm->obsProb = (float**) (buffer + 3);
        status = 3;

    }
    else
    {
        /* check current size */
        int* total= (int*)(((int*)(hmm->obsProb)) - 3);
        int need_size = ( obs_info->obs_x /* * obs_info->obs_y*/ * total_states * sizeof(float) /* +
                           obs_info->obs_y * hmm->num_states * sizeof( CvMatr32f(float*)  )*/ );

        assert( sizeof(float*) == sizeof(int) );

        if ( need_size > (*total) )
        {
            int* buffer = ((int*)(hmm->obsProb)) - 3;
            icvFree( &buffer);
            buffer = (int*)icvAlloc( need_size + 3);
            buffer[0] = need_size;
            buffer[1] = obs_info->obs_y;
            buffer[2] = obs_info->obs_x;

            hmm->obsProb = (float**)(buffer + 3);

            status = 3;
        }
    }
    if (!status)
    {
        int* obsx = ((int*)(hmm->obsProb)) - 1;
        //int* obsy = ((int*)(hmm->obsProb)) - 2;

        assert( /*(*obsy > 0) &&*/ (*obsx > 0) );

        /* is good distribution? */
        if ( (obs_info->obs_x > (*obsx) ) /* || (obs_info->obs_y > (*obsy) ) */ )
            status = 3;
    }

    assert( (status == 0) || (status == 3) );
    /* if bad status - do reallocation actions */
    if ( status )
    {
        float** tmp = hmm->obsProb;
        //float*  tmpf;

        /* distribute pointers of ehmm->obsProb */
/*        for( i = 0; i < hmm->num_states; i++ )
        {
            hmm->u.ehmm[i].obsProb = tmp;
            tmp += obs_info->obs_y;
        }
*/
        //tmpf = (float*)tmp;

        /* distribute pointers of ehmm->obsProb[j] */
/*      for( i = 0; i < hmm->num_states; i++ )
        {
            CvEHMM* ehmm = &( hmm->u.ehmm[i] );

            for( j = 0; j < obs_info->obs_y; j++ )
            {
                ehmm->obsProb[j] = tmpf;
                tmpf += ehmm->num_states * obs_info->obs_x;
            }
        }
*/
        hmm->obsProb = tmp;

    }/* end of pointer distribution */

#if 1
    {
#define MAX_BUF_SIZE  1200
        float  local_log_mix_prob[MAX_BUF_SIZE];
        double local_mix_prob[MAX_BUF_SIZE];
        int    vect_size = obs_info->obs_size;
        CvStatus res = CV_NO_ERR;

        float*  log_mix_prob = local_log_mix_prob;
        double* mix_prob = local_mix_prob;

        int  max_size = 0;
        int  obs_x = obs_info->obs_x;

        /* calculate temporary buffer size */
        //for( i = 0; i < hmm->num_states; i++ )
        //{
        //    CvEHMM* ehmm = &(hmm->u.ehmm[i]);
            CvEHMMState* state = hmm->u.state;

            int max_mix = 0;
            for( j = 0; j < hmm->num_states; j++ )
            {
                int t = state[j].num_mix;
                if( max_mix < t ) max_mix = t;
            }
            max_mix *= hmm->num_states;
            /*if( max_size < max_mix )*/ max_size = max_mix;
        //}

        max_size *= obs_x * vect_size;

        /* allocate buffer */
        if( max_size > MAX_BUF_SIZE )
        {
            log_mix_prob = (float*)icvAlloc( max_size*(sizeof(float) + sizeof(double)));
            if( !log_mix_prob ) return CV_OUTOFMEM_ERR;
            mix_prob = (double*)(log_mix_prob + max_size);
        }

        memset( log_mix_prob, 0, max_size*sizeof(float));

        /*****************computing probabilities***********************/

        /* loop through external states */
        //for( i = 0; i < hmm->num_states; i++ )
        {
        //    CvEHMM* ehmm = &(hmm->u.ehmm[i]);
            CvEHMMState* state = hmm->u.state;

            int max_mix = 0;
            int n_states = hmm->num_states;

            /* determine maximal number of mixtures (again) */
            for( j = 0; j < hmm->num_states; j++ )
            {
                int t = state[j].num_mix;
                if( max_mix < t ) max_mix = t;
            }

            /* loop through rows of the observation matrix */
            //for( j = 0; j < obs_info->obs_y; j++ )
            {
                int  m, n;

                float* obs = obs_info->obs;/* + j * obs_x * vect_size; */
                float* log_mp = max_mix > 1 ? log_mix_prob : (float*)(hmm->obsProb);
                double* mp = mix_prob;

                /* several passes are done below */

                /* 1. calculate logarithms of probabilities for each mixture */

                /* loop through mixtures */
    /*  !!!! */     for( m = 0; m < max_mix; m++ )
                {
                    /* set pointer to first observation in the line */
                    float* vect = obs;

                    /* cycles through obseravtions in the line */
                    for( n = 0; n < obs_x; n++, vect += vect_size, log_mp += n_states )
                    {
                        int k, l;
                        for( l = 0; l < n_states; l++ )
                        {
                            if( state[l].num_mix > m )
                            {
                                float* mu = state[l].mu + m*vect_size;
                                float* inv_var = state[l].inv_var + m*vect_size;
                                double prob = -state[l].log_var_val[m];
                                for( k = 0; k < vect_size; k++ )
                                {
                                    double t = (vect[k] - mu[k])*inv_var[k];
                                    prob -= t*t;
                                }
                                log_mp[l] = MAX( (float)prob, -500 );
                            }
                        }
                    }
                }

                /* skip the rest if there is a single mixture */
                if( max_mix != 1 )
                {
                    /* 2. calculate exponent of log_mix_prob
                          (i.e. probability for each mixture) */
                    res = icvbExp_32f64f( log_mix_prob, mix_prob,
                                            max_mix * obs_x * n_states );
                    if( res < 0 ) goto processing_exit;

                    /* 3. sum all mixtures with weights */
                    /* 3a. first mixture - simply scale by weight */
                    for( n = 0; n < obs_x; n++, mp += n_states )
                    {
                        int l;
                        for( l = 0; l < n_states; l++ )
                        {
                            mp[l] *= state[l].weight[0];
                        }
                    }

                    /* 3b. add other mixtures */
                    for( m = 1; m < max_mix; m++ )
                    {
                        int ofs = -m*obs_x*n_states;
                        for( n = 0; n < obs_x; n++, mp += n_states )
                        {
                            int l;
                            for( l = 0; l < n_states; l++ )
                            {
                                if( m < state[l].num_mix )
                                {
                                    mp[l + ofs] += mp[l] * state[l].weight[m];
                                }
                            }
                        }
                    }

                    /* 4. Put logarithms of summary probabilities to the destination matrix */
                    res = icvbLog_64f32f( mix_prob, (float*)(hmm->obsProb),//[j],
                                            obs_x * n_states );
                    if( res < 0 ) goto processing_exit;
                }
            }
        }

processing_exit:

        if( log_mix_prob != local_log_mix_prob ) icvFree( &log_mix_prob );
        return res;
#undef MAX_BUF_SIZE
    }
#else
/*    for( i = 0; i < hmm->num_states; i++ )
    {
        CvEHMM* ehmm = &(hmm->u.ehmm[i]);
        CvEHMMState* state = ehmm->u.state;

        for( j = 0; j < obs_info->obs_y; j++ )
        {
            int k,m;

            int obs_index = j * obs_info->obs_x;

            float* B = ehmm->obsProb[j];

            // cycles through obs and states
            for( k = 0; k < obs_info->obs_x; k++ )
            {
                CvVect32f vect = (obs_info->obs) + (obs_index + k) * vect_size;

                float* matr_line = B + k * ehmm->num_states;

                for( m = 0; m < ehmm->num_states; m++ )
                {
                    matr_line[m] = icvComputeGaussMixture( vect, state[m].mu, state[m].inv_var,
                                                             state[m].log_var_val, vect_size, state[m].weight,
                                                             state[m].num_mix );
                }
            }
        }
    }
*/
#endif
}


/*F///////////////////////////////////////////////////////////////////////////////////////
//    Name: EstimateTransProb
//    Purpose: The function calculates the state and super state transition probabilities
//             of the model given the images,
//             the state segmentation and the input parameters
//    Context:
//    Parameters: obs_info_array - array of pointers to image observations
//                num_img - length of above array
//                hmm - pointer to HMM structure
//    Returns: void
//
//    Notes:
//F*/
CvStatus icvEstimate1DTransProb( Cv1DObsInfo** obs_info_array,
                                 int num_seq,
                                 CvEHMM* hmm )
{
    int    i, j, k;

    /* as a counter we will use transP matrix */

    /* initialization */

    /* clear transP */
    icvSetZero_32f( hmm->transP, hmm->num_states, hmm->num_states );


    /* compute the counters */
    for (i = 0; i < num_seq; i++)
    {
        int counter = 0;
        Cv1DObsInfo* info = obs_info_array[i];

        for (k = 0; k < info->obs_x; k++, counter++)
        {
            /* compute how many transitions from state to state
               occured */
            int state;
            int nextstate;

            state = info->state[counter];

            if (k < info->obs_x - 1)
            {
                int transP_size = hmm->num_states;

                nextstate = info->state[counter+1];
                hmm->transP[ state * transP_size + nextstate] += 1;
            }
        }
    }

    /* estimate superstate matrix */
    for( i = 0; i < hmm->num_states; i++)
    {
        float total = 0;
        float inv_total;
        for( j = 0; j < hmm->num_states; j++)
        {
            total += hmm->transP[i * hmm->num_states + j];
        }
        //assert( total );

        inv_total = total ? 1.f/total : 0;

        for( j = 0; j < hmm->num_states; j++)
        {
            hmm->transP[i * hmm->num_states + j] =
                hmm->transP[i * hmm->num_states + j] ?
                (float)log( hmm->transP[i * hmm->num_states + j] * inv_total ) : -BIG_FLT;
        }
    }

    return CV_NO_ERR;
}


/*F///////////////////////////////////////////////////////////////////////////////////////
//    Name: MixSegmL2
//    Purpose: The function implements the mixture segmentation of the states of the embedded HMM
//    Context: used with the Viterbi training of the embedded HMM
//
//    Parameters:
//             obs_info_array
//             num_img
//             hmm
//    Returns: void
//
//    Notes:
//F*/
CvStatus icv1DMixSegmL2(CvImgObsInfo** obs_info_array, int num_img, CvEHMM* hmm )
{
    int     k, i, m;

    CvEHMMState* state = hmm->u.state;

    for (k = 0; k < num_img; k++)
    {
        //int counter = 0;
        CvImgObsInfo* info = obs_info_array[k];

        for (i = 0; i < info->obs_x; i++)
        {
            int e_state = info->state[i];
            float min_dist;

            min_dist = icvSquareDistance((info->obs) + (i * info->obs_size),
                                               state[e_state].mu, info->obs_size);
            info->mix[i] = 0;

            for (m = 1; m < state[e_state].num_mix; m++)
            {
                float dist=icvSquareDistance( (info->obs) + (i * info->obs_size),
                                               state[e_state].mu + m * info->obs_size,
                                               info->obs_size);
                if (dist < min_dist)
                {
                    min_dist = dist;
                    /* assign mixture with smallest distance */
                    info->mix[i] = m;
                }
            }
        }
    }
    return CV_NO_ERR;
}

/*F///////////////////////////////////////////////////////////////////////////////////////
//    Name: icvEViterbi
//    Purpose: The function calculates the embedded Viterbi algorithm
//             for 1 image
//    Context:
//    Parameters:
//             obs_info - observations
//             hmm      - HMM
//
//    Returns: the Embedded Viterbi probability (float)
//             and do state segmentation of observations
//
//    Notes:
//F*/
float icvViterbi(Cv1DObsInfo* obs_info, CvEHMM* hmm)
{
    int    i, counter;
    float  log_likelihood;

    //CvEHMMState* first_state = hmm->u.state;

    /* memory allocation for superB */
    /*CvMatr32f superB = picvCreateMatrix_32f(hmm->num_states, obs_info->obs_x );*/

    /* memory allocation for q */
    int* super_q = (int*)icvAlloc( obs_info->obs_x * sizeof(int) );

    /* perform Viterbi segmentation (process 1D HMM) */
    icvViterbiSegmentation( hmm->num_states, obs_info->obs_x,
                            hmm->transP, (float*)(hmm->obsProb), 0,
                            _CV_LAST_STATE, &super_q, obs_info->obs_x,
                             obs_info->obs_x, &log_likelihood );

    log_likelihood /= obs_info->obs_x ;

    counter = 0;
    /* assign new state to observation vectors */
    for (i = 0; i < obs_info->obs_x; i++)
    {
         int state = super_q[i];
         obs_info->state[i] = state;
    }

    /* memory deallocation for superB */
    /*picvDeleteMatrix( superB );*/
    icvFree( &super_q );

    return log_likelihood;
}

CvStatus icvEstimate1DHMMStateParams(CvImgObsInfo** obs_info_array, int num_img, CvEHMM* hmm)

{
    /* compute gamma, weights, means, vars */
    int k, i, j, m;
    int counter = 0;
    int total = 0;
    int vect_len = obs_info_array[0]->obs_size;

    float start_log_var_val = LN2PI * vect_len;

    CvVect32f tmp_vect = icvCreateVector_32f( vect_len );

    CvEHMMState* first_state = hmm->u.state;

    assert( sizeof(float) == sizeof(int) );

    total+= hmm->num_states;

    /***************Gamma***********************/
    /* initialize gamma */
    for( i = 0; i < total; i++ )
    {
        for (m = 0; m < first_state[i].num_mix; m++)
        {
            ((int*)(first_state[i].weight))[m] = 0;
        }
    }

    /* maybe gamma must be computed in mixsegm process ?? */

    /* compute gamma */
    counter = 0;
    for (k = 0; k < num_img; k++)
    {
        CvImgObsInfo* info = obs_info_array[k];
        int num_obs = info->obs_y * info->obs_x;

        for (i = 0; i < num_obs; i++)
        {
            int state, mixture;
            state = info->state[i];
            mixture = info->mix[i];
            /* computes gamma - number of observations corresponding
               to every mixture of every state */
            ((int*)(first_state[state].weight))[mixture] += 1;
        }
    }
    /***************Mean and Var***********************/
    /* compute means and variances of every item */
    /* initially variance placed to inv_var */
    /* zero mean and variance */
    for (i = 0; i < total; i++)
    {
        memset( (void*)first_state[i].mu, 0, first_state[i].num_mix * vect_len *
                                                                         sizeof(float) );
        memset( (void*)first_state[i].inv_var, 0, first_state[i].num_mix * vect_len *
                                                                         sizeof(float) );
    }

    /* compute sums */
    for (i = 0; i < num_img; i++)
    {
        CvImgObsInfo* info = obs_info_array[i];
        int total_obs = info->obs_x;// * info->obs_y;

        float* vector = info->obs;

        for (j = 0; j < total_obs; j++, vector+=vect_len )
        {
            int state = info->state[j];
            int mixture = info->mix[j];

            CvVect32f mean  = first_state[state].mu + mixture * vect_len;
            CvVect32f mean2 = first_state[state].inv_var + mixture * vect_len;

            icvAddVector_32f( mean, vector, mean, vect_len );
            icvAddSquare_32f_C1IR( vector, vect_len * sizeof(float),
                                    mean2, vect_len * sizeof(float), cvSize(vect_len, 1) );
        }
    }

    /*compute the means and variances */
    /* assume gamma already computed */
    counter = 0;
    for (i = 0; i < total; i++)
    {
        CvEHMMState* state = &(first_state[i]);

        for (m = 0; m < state->num_mix; m++)
        {
            int k;
            CvVect32f mu  = state->mu + m * vect_len;
            CvVect32f invar = state->inv_var + m * vect_len;

            if ( ((int*)state->weight)[m] > 1)
            {
                float inv_gamma = 1.f/((int*)(state->weight))[m];

                icvScaleVector_32f( mu, mu, vect_len, inv_gamma);
                icvScaleVector_32f( invar, invar, vect_len, inv_gamma);
            }

            icvMulVectors_32f(mu, mu, tmp_vect, vect_len);
            icvSubVector_32f( invar, tmp_vect, invar, vect_len);

            /* low bound of variance - 0.01 (Ara's experimental result) */
            for( k = 0; k < vect_len; k++ )
            {
                invar[k] = (invar[k] > 0.01f) ? invar[k] : 0.01f;
            }

            /* compute log_var */
            state->log_var_val[m] = start_log_var_val;
            for( k = 0; k < vect_len; k++ )
            {
                state->log_var_val[m] += (float)log( invar[k] );
            }

            state->log_var_val[m] *= 0.5;

            /* compute inv_var = 1/sqrt(2*variance) */
            icvScaleVector_32f(invar, invar, vect_len, 2.f );
            icvbInvSqrt_32f(invar, invar, vect_len );
        }
    }

    /***************Weights***********************/
    /* normilize gammas - i.e. compute mixture weights */

    //compute weights
    for (i = 0; i < total; i++)
    {
        int gamma_total = 0;
        float norm;

        for (m = 0; m < first_state[i].num_mix; m++)
        {
            gamma_total += ((int*)(first_state[i].weight))[m];
        }

        norm = gamma_total ? (1.f/(float)gamma_total) : 0.f;

        for (m = 0; m < first_state[i].num_mix; m++)
        {
            first_state[i].weight[m] = ((int*)(first_state[i].weight))[m] * norm;
        }
    }

    icvDeleteVector( tmp_vect);
    return CV_NO_ERR;
}


#endif