sammon.cpp

face recognition test source code

/************************************************************************
 *                                                                      *
 *  Program packages 'lvq_pak' and 'som_pak' :                          *
 *                                                                      *
 *  sammon.c                                                            *
 *  - generates a Sammon mapping from a given list                      *
 *                                                                      *
 *  Version 3.0                                                         *
 *  Date: 1 Mar 1995                                                    *
 *                                                                      *
 *  NOTE: This program package is copyrighted in the sense that it      *
 *  may be used for scientific purposes. The package as a whole, or     *
 *  parts thereof, cannot be included or used in any commercial         *
 *  application without written permission granted by its producents.   *
 *  No programs contained in this package may be copied for commercial  *
 *  distribution.                                                       *
 *                                                                      *
 *  All comments  concerning this program package may be sent to the    *
 *  e-mail address '[e-mail]'.                                *
 *                                                                      *
 ************************************************************************/

#include "header.h"
#include <stdio.h>
#include <math.h>
#include <float.h>
#include <string.h>
#include "matlvq.h"
#include "lvq.h"
#include "lvq_pak.h"
#include "tools.h"
#include "labels.h"
#include "sammon.h"
  
/*#define TRUE  1
#define FALSE 0*/
#define MAGIC 0.2
  
#define max(x,y) ((x) > (y) ? (x) : (y))
#define min(x,y) ((x) < (y) ? (x) : (y))

struct entries *sammon_iterate(struct entries *codes, double length, float *error)
{
  int i, j, k;
  long noc = 0;
  float e1x, e1y, e2x, e2y;
  float dpj;
  float dq, dr, dt;
  struct data_entry *entr, *entr1, *prev, tmp;
  float *x, *y;
  float *xu, *yu, *dd;
  float xd, yd;
  float xx, yy;
  float e, tot;
  int mutual;
  float d, ee;
  DIST_FUNCTION *distance = vector_dist_euc; /* tama jotenkin parametrina */
  eptr p1, p2;
  struct entries *newent;
  int dim = codes->dimension;

  /* How many entries? */
  noc = codes->num_entries;

  /* Allocate dynamical memory */
  x = (float *) oalloc(sizeof(float) * noc);
  y = (float *) oalloc(sizeof(float) * noc);
  xu = (float *) oalloc(sizeof(float) * noc);
  yu = (float *) oalloc(sizeof(float) * noc);
  dd = (float *) oalloc(sizeof(float) * (noc * (noc - 1) / 2));

  /* Initialize the tables */
  for (i = 0; i < noc; i++) {
    x[i] = (float) (orand() % noc) / noc;
    y[i] = (float) (i) / noc;
  }

  /* Compute the mutual distances between entries */
  mutual = 0;
  entr = rewind_entries(codes, &p1);
  entr = next_entry(&p1);
  while (entr != NULL) {
    entr1 = rewind_entries(codes, &p2);
    while (entr1 != entr) {
      dd[mutual] = distance(entr, entr1, dim);
      /*if (dd[mutual] == 0.0) {
	fprintf(stderr, "Identical entries in codebook\n");
      } */
      mutual++;
      entr1 = next_entry(&p2);
    }
    entr = next_entry(&p1);
  }
  
  /* Iterate */              
  e = length+1;
  for (i = 0; 1/*i < length*/; i++) {
  		if (length > 1.0 && i > length)
  			break;
  		if (length <= 1.0 && e < length)
  			break;
    for (j = 0; j < noc; j++) {
      e1x = e1y = e2x = e2y = 0.0;
      for (k = 0; k < noc; k++) {
	if (j == k)
	  continue;
	xd = x[j] - x[k];
	yd = y[j] - y[k];
	dpj = (float) sqrt((double) xd * xd + yd * yd);

	/* Calculate derivatives */
	if (k > j)
	  dt = dd[k * (k - 1) / 2 + j];
	else
	  dt = dd[j * (j - 1) / 2 + k];
	
	dq = dt - dpj;
	dr = dt * dpj;
	e1x += xd * dq / dr;
	e1y += yd * dq / dr;
	e2x += (dq - xd * xd * (1.0 + dq / dpj) / dpj) / dr;
	e2y += (dq - yd * yd * (1.0 + dq / dpj) / dpj) / dr;
      }
      /* Correction */
      xu[j] = x[j] + MAGIC * e1x / fabs(e2x);
      yu[j] = y[j] + MAGIC * e1y / fabs(e2y);
    }
    
    /* Move the center of mass to the center of picture */
    xx = yy = 0.0;
    for (j = 0; j < noc; j ++) { 
      xx += xu[j];
      yy += yu[j];
    }
    xx /= noc;
    yy /= noc;
    for (j = 0; j < noc; j ++) {
      x[j] = xu[j] - xx;
      y[j] = yu[j] - yy;
    }
    
    /* Error in distances */
    e = tot = 0.0;
    mutual = 0;
    for (j = 1; j < noc; j ++)
      for (k = 0; k < j; k ++) {
	d = dd[mutual];
	tot += d;
	xd = x[j] - x[k];
	yd = y[j] - y[k];
	ee = d - (float) sqrt((double) xd * xd + yd * yd);
	e += (ee * ee / d);
	mutual++;
      }
    if (tot != 0.0)
    	e /= tot;
    ifverbose(2)
      Msg(0, "Mapping error: %7.3f\n", e);
    /*if (verbose(-1) == 1)
      mprint((long) (length-i));*/
  }
  /*if (verbose(-1) == 1)
    mprint((long) 0);
  if (verbose(-1) == 1)
    fprintf(stderr, "\n");*/

  *error = e;
  
  newent = alloc_entries();
  newent->dimension = 2;
  newent->xdim = codes->xdim;
  newent->ydim = codes->ydim;
  newent->topol = codes->topol;
  newent->neigh = codes->neigh;

  /* Copy the data to return variable */

  prev = &tmp;
  prev->next = NULL;
  for (i = 0, entr = rewind_entries(codes, &p1); i < noc; i++, entr = next_entry(&p1)) {
    entr1 = init_entry(newent, NULL);
    prev->next = entr1;
    prev = entr1;
    entr1->points[0] = x[i];
    entr1->points[1] = y[i];
    copy_entry_labels(entr1, entr);
  }
  newent->entry = tmp.next;
  newent->num_entries = noc;

  return(newent);
}

struct entries *remove_identicals(struct entries *codes, int *rem)
{
  struct data_entry *entr, *entr1, *prev;
  DIST_FUNCTION *distance = vector_dist_euc; /* tama parameettereista */
  eptr p1, p2;
  int dim = codes->dimension;
  int ii, ij;
  
  /* Compute the mutual distances between entries */
  /* Remove the identical entries from the list */
  *rem = 0;
  entr = rewind_entries(codes, &p1);
  ii = 1; ij = 1;
  p2.parent = p1.parent;
  while (entr != NULL) {
    p2.current = p1.current;
    p2.index = p1.index;
    entr1 = next_entry(&p2);
    ij = ii + 1;
    prev = p1.current;
    while (entr1 != NULL) {
      if (distance(entr, entr1, dim) == 0.0) {
	/*fprintf(stderr, "Identical entries in codebook ");
        fprintf(stderr, "(entries %d, %d), removing one.\n", ii, ij);*/
        *rem = 1;
	codes->num_entries--;
	prev->next = entr1->next;
	p2.current = prev;
	free_entry(entr1);
	entr1 = next_entry(&p2);
        ij++; ij++;
      }
      else {
	prev = entr1;
        entr1 = next_entry(&p2);
        ij++;
      }
    }
    entr = next_entry(&p1);
    ii++;
  }
  
  return(codes);
}

int put_sammon(struct teach_params *tparams, char *name, long num, int total)
	{
	int ret;
	float error;
  	struct entries *newent;
	CMatrix *pmat = NULL;
	struct entries *newcodes;
	int rem;
	struct data_entry *entr, *new_entr, *prev;
	eptr p;
	PSAMMON_PICT pict = new SAMMON_PICT;
	
	if (tparams->sammon == NULL)
		return -1;
	
	/*copy the codebook vectors into a new place in order to be able to remove the identical*/
	/*codebook vectors without having to cause trouble in the rest of the process*/
   newcodes = copy_entries(tparams->codes);
   if (newcodes == NULL)
   	return -1;
	entr = rewind_entries(tparams->codes, &p);
	prev = NULL;
	while (entr != NULL) 
		{
		new_entr = copy_entry(newcodes, entr); 
		if (new_entr == NULL) 
			return -1;
		if (prev != NULL)
			prev->next = new_entr;
		else
			{
			newcodes->current = new_entr;
			newcodes->entry = new_entr;
			}
		prev = new_entr;
		new_entr->next = NULL;
		entr = next_entry(&p);
      newcodes->num_entries++;
      newcodes->num_loaded++;
      }

   remove_identicals(newcodes, &rem);

	newent = sammon_iterate(newcodes, tparams->sammon->length, &error);
	if (newent == NULL)
		{
		close_entries(newcodes);
		return -1;
		}
	if (error > 1.0)
		{	/*Mapping error, drop the picture*/
		Msg(0, "Mapping error (error = %f), sammon view of the learning it. num %.0f of %s dropped",
			error, (double)num, name);
		close_entries(newcodes);
		close_entries(newent);
		return 0;
		}
	
	/*free the new codebook vectors*/
	close_entries(newcodes);
		
	ret = EntriestoMat(&pmat, newent);
	if (ret < 0)
		{
		close_entries(newent);
		return -1;
		}
	
 	pict->pmat = pmat;
 	strcpy(pict->name, name);
 	pict->lnum = num;
 	pict->total = total;
 	pict->error = error;
 	tparams->sammon->AddPict(pict);
	
	close_entries(newent);
	return 0;
	}