sparse.c

射线追踪程序

/*
 * This file is part of tomo3d
 *
 * Copyright (C) 2002, 2003, Sebastien Judenherc <[e-mail]>
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307
 * USA
 *
 */
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <assert.h>
#include <math.h>
#include "sparse.h"

#define EPS 1.0e-6

#define OK fprintf(stderr,"FILE %s LINE %d\n",__FILE__,__LINE__);

int writematrix(f,tmp)
FILE *f;
struct matrix_t * tmp;
{ float *line;
  int i,j;

  fprintf(stderr,"Writing matrix ...");
  fwrite(&(tmp->nl),sizeof(int),1,f);
  fwrite (&(tmp->nc),sizeof(int),1,f);
  fprintf(stderr,"%d:%d...",tmp->nl,tmp->nc);
  line=calloc(tmp->nc,sizeof(float));
  if (!line) {
     free(tmp);
     return(2);
  }
  for (i=0; i<tmp->nl; i++) {
      for (j=0; j<tmp->nc; j++) sparsegetval(tmp, i, j, &(line[j]));
      fwrite(line,sizeof(float),tmp->nc,f);
  }
  free(line);
  fprintf(stderr,"OK\n");
  return(0);
}

struct matrix_t * readmatrix(f)
FILE *f;
{ struct matrix_t * tmp;
  float *line;
  int L,C,i,j;

  fprintf(stderr,"Reading matrix ...");
  fread(&L,sizeof(int),1,f);
  fread(&C,sizeof(int),1,f);
  fprintf(stderr,"%d:%d...",L,C);
  tmp=newmatrix(L,C);
  if (!tmp) {
     return(NULL);
  }
  line=calloc(C,sizeof(float));
  if (!line) {
     free(tmp);
     return(NULL);
  }
  for (i=0; i<L; i++) {
      fread(line,sizeof(float),C,f);
      for (j=0; j<C; j++)
          if (fabs(line[j])>EPS)
             sparseputval(tmp, i, j, line[j]);
  }
  free(line);
  fprintf(stderr,"OK\n");
  return(tmp);
}

struct matrix_t * newmatrix(nl,nc)
int nl, nc;
{ struct matrix_t * tmp;
  int i;
  tmp=(struct matrix_t *)calloc(1,sizeof(struct matrix_t));
  if (!tmp) return(NULL);

  tmp->nl=nl;
  tmp->nc=nc;
  tmp->M=(struct line_t *)calloc(nl,sizeof(struct line_t));
  if (!tmp->M) {
     free(tmp);
     return(NULL);
  }
  for (i=0; i<nl; i++) tmp->M[i].index=-1;
  return(tmp);
}

void sparsefreematrix(mat)
struct matrix_t * mat;
{ int i;
  struct line_t * cur, *last, *tmp;
  for (i=0; i<mat->nl; i++) {
      tmp=mat->M[i].next;
      while (tmp->next) {
            for (last=tmp,cur=tmp->next; cur->next; last=cur,cur=cur->next) ;
            last->next=NULL;
            if (cur) free(cur);
            cur=NULL;
      }
      if (tmp->next) free(tmp->next);
      tmp->next=NULL;
  }
  free(mat->M);
  free(mat);
}

int sparsedelval(mat, i, j)
struct matrix_t * mat;
int i,j;
{ struct line_t *cur,*tmp;

  cur=tmp=&(mat->M[i]);

  /* la ligne est elle vide ? */
  if (cur->index==-1) return(1); /* coup pour rien */

  /* ou sa place est elle la premiere ? */
  if (cur->index==j) {
     if (cur->next) { /* mais pas la derniere */
        tmp=cur->next;
        memcpy(cur,tmp,sizeof(struct line_t));
        free(tmp);
        return(0);
     } else { cur->index=-1; /* la ligne devient vide */
              cur->next=NULL;
              return(0);
     }
  }

  /* on cherche sa place */
  for (tmp=cur,cur=cur->next; cur; tmp=cur,cur=cur->next)
      if (j>=cur->index) break;

  /* il n'existe pas: ok */
  if (!cur) return(0);

  /* il existe, on le tue */
  if (j==cur->index) {
     tmp->next=cur->next;
     free(cur);
     return(0);
  }

  return(1);
}


int sparseputval(mat, i, j, val)
struct matrix_t * mat;
int i,j;
float val;
{ struct line_t *cur,*tmp, *tmp2;

  if (fabs(val)<EPS)
     return(sparsedelval(mat, i, j));

  cur=tmp=&(mat->M[i]);

  /* est-ce le premier et seul non nul de la ligne ? */
  /* ou sa place est elle la premiere ? */
  if ((cur->index==-1)||(cur->index==j)) {
     cur->val=val;
     cur->index=j;
     return(1);
  }

  /* on cherche sa place */
  for (; cur; tmp=cur,cur=cur->next) if (j>cur->index) break;

  if (cur) {
          /* est-ce le premier de la ligne ? */
     if (cur==&(mat->M[i])) {
        tmp=(struct line_t *)calloc(1,sizeof(struct line_t));
        if (!tmp) return(-1);
        memcpy(tmp,cur,sizeof(struct line_t));
   
        cur->index=j;
        cur->val=val;
        cur->next=tmp;
        return(0);
     }

     /* la bonne place est entre tmp et cur (cur=tmp->next) */
     tmp2=(struct line_t *)calloc(1,sizeof(struct line_t));
     tmp2->index=j;
     tmp2->val=val;
     tmp2->next=cur;
     tmp->next=tmp2;
     return(0);
  }
       /* est-ce le dernier de la liste ? */
  if (cur==NULL) {
     tmp->next=(struct line_t *)calloc(1,sizeof(struct line_t));
     cur=tmp->next;
     if (!cur) return(-1);
     cur->index=j;
     cur->val=val;
     return(0);
  }
  return(-1); /* ??????? */
}

int sparsegetval( mat, i, j, val )
struct matrix_t * mat;
int i, j;
float *val;
{ struct line_t *cur;
  *val=0.0;

  cur=&(mat->M[i]);
  if (cur->index==-1) return(0);

  for ( ; cur; cur=cur->next)
      if (j>=cur->index) break;

  if (cur) {
          /* on a trouve sa place exacte ? */
     if (cur->index==j) {
        *val=cur->val;
        return(1);
     }
  }
  return(0);
}