#include <proeikon>
#include <imgpv>
#include <vector>
#include <queue>

bool ignora255=true;
// Ignora pixels de q com valor 255 - usado para invariancia afim
bool linhafina=true;

void Ignora255(IMGGRY qu, IMGGRY& valido, int tnl)
// So e' usado para Aforapro
{ qu.backg()=255; //IMGGRY qushow=qu; 
  int raio=tnl/2;
  if (qu.nl()-tnl+1!=valido.nl() || qu.nc()-tnl+1!=valido.nc()) erro("Erro Ignora255: dimensoes incompativeis");
  for (int l=0; l<valido.nl(); l++)
    for (int c=0; c<valido.nc(); c++) {
      for (int l2=-raio; l2<=raio; l2++)
        for (int c2=-raio; c2<=raio; c2++) {
          double r=sqrt(l2*l2+c2*c2);
          if (int(r)<=raio && qu(l+l2+raio,c+c2+raio)==255) { 
            valido(l,c)=255; 
            //qushow(l+raio,c+raio)=0; 
          }
        }
    }
  //pvMostra(qushow);
}

void impmat(MATRIZ<PONTOI2> a)
// MATRIZ e' matriz.
// PONTOI2 e' um vetor de duas dimensoes, cada elemento e' inteiro.
// PONTOI2 p; p.l()=3; p[0]=3; p.x()=3;
{ for (int l=0; l<a.nl(); l++) {
    for (int c=0; c<a.nc(); c++)
      printf("(%d,%d) ",a(l,c).l(),a(l,c).c());
    printf("\n");
  }
}

//const double minampl=0.01;
const double minampl=0.0025; // ????
struct FEA { // features
  double rotmap; // - 1 = rotmap: Numeros entre -pi a +pi. orientacao canonica
  FUNCAO rotvet; // angulos - nrad-1 = rotvet: vetor de angulos entre 0 e +2pi
                 // vector of radial angles
  FUNCAO radvet; // comprimentos - nrad - radvet: Vetor de comprimento 1
                 // vector of radial magnitudes
  FUNCAO cirvet; // 2*ncir - cirvet: Vetor de comprimento 1
                 // vector of circular features
  // FUNCAO = VETOR<double>
};
typedef MATRIZ<FEA> MFEA;

struct PESOTOTAL {
  double vra_wt,wm,wa,wc,wt;
  PESOTOTAL(int nrad, int ncir) {
    vra_wt=0.0;
    for (int i=0; i<nrad-1; i++) // kth=i+2
      vra_wt += 1.0/(i+2.0);
    wm=nrad-1;
    wa=nrad-1;
    wc=2*ncir-1;
    wt=wm+wa+wc;
    wm=wm/wt;
    wa=wa/wt;
    wc=wc/wt;
  } // equacao (20) do artigo
};

inline double rpeso(double r, int raio)
// peso de nucleo radial
// primeiro fator da equacao (7)
{ double fim=raio+0.5;
  double inicio=1.0;
  if (r<=inicio || fim<=r) return 0.0;
  else return sqrt((r-inicio)*(fim-r));
}

IMGCPX gerarqux(int kth, int qnl, bool print_qux=false)
// Gera um nucleo radial de ordem kth e diametro qnl.
// Cada elemento de IMGCPX e' complex<float>
// complex<float> == CPX
// complex<double> == CPXD
{ static CPXD jx(0.0, 1.0);
  IMGCPX qux(qnl,qnl);
  int raio=qnl/2;

  double total=0.0;
  for (int l=qux.minl(); l<=qux.maxl(); l++)
    for (int c=qux.minc(); c<=qux.maxc(); c++) {
      double r=sqrt(l*l+c*c);
      double peso=rpeso(r,raio);
      total += peso;
      if (peso>0.0) {
        //double angulo=-kth*arg(CPXD(c,-l));
        double angulo=kth*arg(CPXD(c,-l));
        qux.atc(l,c) = CPX(peso*cos(angulo),peso*sin(angulo));
      } else {
        qux.atc(l,c) = CPX(0.0,0.0);
      }
    }
  for (int l=qux.minl(); l<=qux.maxl(); l++)
    for (int c=qux.minc(); c<=qux.maxc(); c++)
      if (qux.atc(l,c)!=CPX(0.0,0.0))
        qux.atc(l,c) /= total;
  if (print_qux) imp(qux,"rqux.txt");

  return qux;
}

IMGCPX geracqux(int kth, int qnl, bool print_qux=false)
// Gera um nucleo circular de ordem kth e diametro qnl.
// Cada elemento de IMGCPX e' complex<float>
{ static CPXD jx(0.0, 1.0);
  IMGCPX qux(qnl,qnl);
  int raio=qnl/2;

  // continuo
  double total=0.0;
  for (int l=qux.minl(); l<=qux.maxl(); l++)
    for (int c=qux.minc(); c<=qux.maxc(); c++) {
      double r=sqrt(l*l+c*c);
      if (int(r)>raio) qux.atc(l,c)=CPX(0.0,0.0);
      else {
        double peso;
        if (r>0.0) peso=1.0/(M_2PI*r);
        else peso=0.73;
        total += peso;
        //double angulo = -kth * (r/(raio+1)) * M_2PI;
        double angulo = kth * (r/(raio+1)) * M_2PI;
        qux.atc(l,c) = CPX(peso*cos(angulo),peso*sin(angulo));
      }
    }
  for (int l=qux.minl(); l<=qux.maxl(); l++)
    for (int c=qux.minc(); c<=qux.maxc(); c++)
      if (qux.atc(l,c)!=CPX(0.0,0.0))
        qux.atc(l,c) /= total;
  if (print_qux) imp(qux,"cqux.txt");
  return qux;
}

void geraqux(I3DCPX& qux, int nrad, int ncir, int tnl)
// I3DCPX e' como se fosse vetor de IMGCPX.
// I3DCPX z: z(s,l,c) z(s)
// Chamar gerarqux nrad vezes. Depois, chama geracqux ncir vezes
// tnl = diametro do template
// Exemplo: nrad=3 e ncir=2
// qux tera 5 slices: 0..2 referem-se a nucleos radiais
//                    3..4 referem-se a nucleos circulares
{ qux.resize(nrad+ncir,tnl,tnl);
  for (int i=0; i<nrad; i++) {
    int ip=i+1;
    qux(i)=gerarqux(ip,tnl);
  }
  for (int i=nrad; i<nrad+ncir; i++) {
    int ip=i-nrad+1;
    qux(i)=geracqux(ip,tnl);
  }
}

void convqux(IMGGRY& qu, I3DCPX qux, I3DCPX& quc)
// Convolucao de qu com cada nucleo de qux, gerando quc
// Se qux tiver 5 slices, quc tambem tera 5 slices
// Cuidado: Numero de linhas e colunas de quc e' menor que os de qu
// Nota: Na nova versao, seria melhor usar convolucao modo "same"
{ int tnl=qux.nl();
  int nl=qu.nl()-tnl+1;
  int nc=qu.nc()-tnl+1;
  quc.resize(qux.ns(),nl,nc);

  IMGCPX qx=qu;

  // Reflexao de qux
  for (int i=0; i<qux.ns(); i++)
    qux(i)=~qux(i);

  pvMultConvV(qx, qux, quc); // modo valid
}

void gerafea(I3DCPX& quc, MFEA& mfea, IMGGRY& valido, int nrad, int ncir, int tnl)
// Gera mfea (matriz de features) a partir do quc
// e gera imagem valido indica pixels estaveis e instaveis
//   valido(l,c)==255 indica nao-valido
// Nota: Na nova versao, teria que achar os pixels "melhores", nao apenas estaveis
{
  mfea.resize(quc.nl(),quc.nc());
  //rotmap.resize(quc.nl(),quc.nc()); // dim 1 (orientacao canonica)
  //rotvet.resize(quc.nl(),quc.nc()); // dim nrad-1 (vector of radial angles)
  //radvet.resize(quc.nl(),quc.nc()); // dim nrad (vector of radial magnitudes)
  //cirvet.resize(quc.nl(),quc.nc()); // dim 2*ncir (vector of circular features)

  valido.resize(quc.nl(),quc.nc()); valido.fill(BYTE(0));
  int n=quc(0).n(); // numero de elementos de quc, mfea e valido

  for (int i=0; i<n; i++) { // rotmap: Numeros entre -pi a +pi
    CPX x=quc(0).atf(i);
    if (abs(x)<minampl) valido.atf(i)=255;
    mfea.atf(i).rotmap=arg(x); // arg = fase do numero complexo
  }

  for (int i=0; i<n; i++) { // rotvet: vetor de angulos entre 0 e +2pi
    FUNCAO& v=mfea.atf(i).rotvet;
    v.resize(nrad-1);
    double a1=mfea.atf(i).rotmap;
    for (int s=1; s<nrad; s++) {
      int kth=s+1;
      CPX x=quc(s).atf(i);
      double dk=abs(x);
      if (dk<minampl) valido.atf(i)=255;
      double ak=arg(x);
      v.atf(s-1)=modulo((ak-kth*a1),M_2PI); // OK
    }
  }

  for (int i=0; i<n; i++) { // radvet: Vetor de comprimento 1
    FUNCAO& v=mfea.atf(i).radvet;
    v.resize(nrad);
    for (int s=0; s<v.n(); s++) v.atf(s)=abs(quc(s).atf(i));
    v=versorl1(v);
  }

  if (ncir>0) {
    for (int i=0; i<n; i++) { // cirvet: Vetor de comprimento 1
      FUNCAO& v=mfea.atf(i).cirvet;
      v.resize(2*ncir);
      for (int s=0; s<ncir; s++) {
        v.atf(2*s)  =real(quc(s+nrad).atf(i));
        v.atf(2*s+1)=imag(quc(s+nrad).atf(i));
      }
      v=versorl1(v);
    }
  }
}

inline double difrotvet(FUNCAO& a, FUNCAO& b, PESOTOTAL& pesototal) // Entre 0 e 1
// kth radarg recebe peso 1/kth: 1/2, 1/3, 1/4, etc.
// Equacao (15) do artigo
{ //double total=0.0;
  //for (int i=0; i<a.n(); i++) // kth=i+2
  //  total += 1.0/(i+2.0);

  double soma=0.0;
  for (int i=0; i<a.n(); i++)
    soma += (1.0/((i+2.0)*pesototal.vra_wt)) * difangrad(a.atf(i),b.atf(i));
  return soma/M_PI;
}

inline double diffea(FEA& f1, FEA& f2, PESOTOTAL& pesototal)
// Nota: Esta funcao esta' reescrito em procura1.
// Se alterar esta funcao, tem que alterar tambem procura1.
// Equacao (20) do artigo
{ double t1=difrotvet(f1.rotvet,f2.rotvet,pesototal); // Entre 0 e 1
  double t2=0.5*distancial1(f1.radvet,f2.radvet);       // Entre 0 e 1
  double t3=0.5*distancial1(f1.cirvet,f2.cirvet);       // Entre 0 e 1

  //if (t1>1.0)  { printf("%f\n",t1); erro("Erro: t1>1.0"); }
  //if (t2>1.0)  { printf("%f\n",t2); erro("Erro: t2>1.0");}
  //if (t3>1.0)  { printf("%f\n",t3); erro("Erro: t3>1.0"); }

  return (pesototal.wa*t1+pesototal.wm*t2+pesototal.wc*t3);
}

void stablecand(MFEA& mfea, IMGGRY& valido, double precisao, PESOTOTAL& pesototal)
// Recebe mfea e valido (0==valido 255==angulo-indeterminado)
// classifica os pixels instaveis como nao-validos
{
  for (int l=0; l<valido.nl(); l++) {
    valido.atf(l,0)=255;
    valido.atf(l,valido.nc()-1)=255;
  }
  for (int c=0; c<valido.nc(); c++) {
    valido.atf(0,c)=255;
    valido.atf(valido.nl()-1,c)=255;
  }

  IMGGRY validoe=dilat(valido,3,3); // Pixels que tem os 8 vizinhos validos.
  valido=validoe;

  for (int l=1; l<validoe.nl()-1; l++) // rotmap
    for (int c=1; c<validoe.nc()-1; c++) {
      if (validoe.atf(l,c)) {
        double maxdif=-infinito;
        for (int l2=-1; l2<=1; l2++)
          for (int c2=-1; c2<=1; c2++) {
            double t=difangrad(mfea.atf(l,c).rotmap,mfea.atf(l+l2,c+c2).rotmap)/M_PI;
            if (t>maxdif) maxdif=t;
          }
        //if (maxdif>(M_2PI/3.6)*precisao) valido.atf(l,c)=255; // Original - nao dividia digangrad
        if (maxdif>precisao) valido.atf(l,c)=255;
      }
    }

  for (int l=1; l<validoe.nl()-1; l++) // fea
    for (int c=1; c<validoe.nc()-1; c++) {
      if (validoe.atf(l,c)) {
        double maxdif=-infinito;
        for (int l2=-1; l2<=1; l2++)
          for (int c2=-1; c2<=1; c2++) {
            double t=diffea(mfea.atf(l,c),mfea.atf(l+l2,c+c2),pesototal);
            if (t>maxdif) maxdif=t;
          }
        if (maxdif>precisao) valido.atf(l,c)=255;
      }
    }
}

IMGFLT distgeo(IMGBIN a)
// Calcula distancia geodesica 4-conexa
// Aparentemente, esta funcao contem erro.
// Verificacao de pixel preto tem que ser feita apos retirar da fila
{ IMGBIN b=-(a^(a+IMGBIN(3,3, 0,1,0, 1,1,1, 0,1,0)));
  queue<int> q;
  for (int l=0; l<a.nl(); l++)
    for (int c=0; c<a.nc(); c++)
      if (b(l,c)==preto) { q.push(l); q.push(c); q.push(1); a(l,c)=branco; }

  IMGFLT g(a.nl(),a.nc(),0.0);
  while (!q.empty()) {
    int l=q.front(); q.pop();
    int c=q.front(); q.pop();
    int k=q.front(); q.pop();
    g(l,c)=k; k++;
    if (a(l-1,c)==preto) { q.push(l-1); q.push(c); q.push(k); a(l-1,c)=branco; }
    if (a(l+1,c)==preto) { q.push(l+1); q.push(c); q.push(k); a(l+1,c)=branco; }
    if (a(l,c-1)==preto) { q.push(l); q.push(c-1); q.push(k); a(l,c-1)=branco; }
    if (a(l,c+1)==preto) { q.push(l); q.push(c+1); q.push(k); a(l,c+1)=branco; }
  }
  return g;
}

void relogioretangular2(IMGCOR& a, int l, int c, double angrad, int nl, int nc, int q=-1, int s=-1)
{ double co=cos(angrad); double si=sin(angrad);
  int lc=(nl-1)/2; int cc=(nc-1)/2;

  int dl0=arredonda(-co*lc             - si*cc);
  int dc0=arredonda( si*lc             - co*cc);
  int dl1=arredonda(-co*lc             + si*(nc-cc-1) );
  int dc1=arredonda( si*lc             + co*(nc-cc-1) );
  int dl2=arredonda( co*(nl-lc-1) - si*cc);
  int dc2=arredonda(-si*(nl-lc-1) - co*cc);
  int dl3=arredonda( co*(nl-lc-1) + si*(nc-cc-1) );
  int dc3=arredonda(-si*(nl-lc-1) + co*(nc-cc-1) );
  int dl=arredonda(-co*lc);
  int dc=arredonda( si*lc);

  if (linhafina==false) {
    thickline(a,l+dl0,c+dc0,l+dl1,c+dc1,COR(0,0,255),3);
    thickline(a,l+dl1,c+dc1,l+dl3,c+dc3,COR(0,0,255),3);
    thickline(a,l+dl3,c+dc3,l+dl2,c+dc2,COR(0,0,255),3);
    thickline(a,l+dl2,c+dc2,l+dl0,c+dc0,COR(0,0,255),3);
    thickline(a,l,c,l+dl,c+dc,COR(0,0,255),3);
    ponto(a,l,c,5,COR(0,255,0));
  } else {
    reta(a,l+dl0,c+dc0,l+dl1,c+dc1,COR(0,0,255));
    reta(a,l+dl1,c+dc1,l+dl3,c+dc3,COR(0,0,255));
    reta(a,l+dl3,c+dc3,l+dl2,c+dc2,COR(0,0,255));
    reta(a,l+dl2,c+dc2,l+dl0,c+dc0,COR(0,0,255));
    reta(a,l,c,l+dl,c+dc,COR(0,0,255));
    ponto(a,l,c,3,COR(0,255,0));
  }

  static char st[256];
  if (q>=0) {
    sprintf(st,"%d",q);
    puttxt(a,l-10,c+2,st,COR(0,255,0));
  }
  if (s>=0) {
    sprintf(st,"%d",s);
    puttxt(a,l+10,c+2,st,COR(0,255,0));
  }
}

void procura(MFEA& anf, FEA& tef, IMGFLT& mat, PESOTOTAL& pesototal)
// Recebe anf e tef. 
// Calcula a diferenca entre anf(l,c) e tef em cada pixel
// Armazena o resultado no mat
{ mat.resize(anf.nl(),anf.nc()); //mat.fill(infinito);
  for (int i=0; i<anf.n(); i++)
    mat.atf(i)=diffea(tef,anf.atf(i),pesototal);
}

void procura1(MFEA& anf, FEA& f1, IMGFLT& mat, double deteccao, PESOTOTAL& pesototal)
// versao acelerada de procura quando se deseja achar um unico matching
{ mat.resize(anf.nl(),anf.nc()); //mat.fill(infinito);

  //double total=0.0;
  //for (int i=0; i<f1.cirvet.n(); i++) // kth=i+2
  //  total += 1.0/(i+2.0);

  double soma,t3,t2,t1;
  for (int i=0; i<anf.n(); i++) {
    FEA& f2=anf.atf(i);
    soma=0.0;

    t2=0.0;
    for (int j=0; j<f1.radvet.n(); j++)
      t2 += abs(f1.radvet.atf(j)-f2.radvet.atf(j));
    soma = pesototal.wm * 0.5 * t2;
    if (soma>deteccao) { mat.atf(i)=1.0; continue; }

    t3=0.0;
    for (int j=0; j<f1.cirvet.n(); j++)
      t3 += abs(f1.cirvet.atf(j)-f2.cirvet.atf(j));
    soma += pesototal.wc * 0.5 * t3;
    if (soma>deteccao) { mat.atf(i)=1.0; continue; }

    t1=0.0;
    for (int j=0; j<f1.rotvet.n(); j++)
      t1 += (1.0/((j+2.0)*pesototal.vra_wt)) * difangrad(f1.rotvet.atf(j),f2.rotvet.atf(j));
    soma += pesototal.wa * t1 / M_PI;
    if (soma>deteccao) { mat.atf(i)=1.0; continue; }

    mat.atf(i)=soma;
  }
}

void imp(MFEA& mfea, string nome, int tnl)
{ FILE* arq=fopen(nome.c_str(),"wb");

  char st[]="IMGFEA*\x1a";
  fwrite(st,1,8,arq);

  int n=mfea.nl();
  fwrite(&n,sizeof(n),1,arq);
  n=mfea.nc();
  fwrite(&n,sizeof(n),1,arq);
  int nrad=mfea(0).radvet.n(); //nrad
  fwrite(&nrad,sizeof(nrad),1,arq);
  int ncir=mfea(0).cirvet.n()/2; //ncir
  fwrite(&ncir,sizeof(ncir),1,arq);
  fwrite(&tnl,sizeof(tnl),1,arq);

  for (int i=0; i<mfea.n(); i++) {
    fwrite(&(mfea(i).rotmap),  sizeof(double), 1,      arq);
    fwrite(mfea(i).rotvet.vet, sizeof(double), nrad-1, arq);
    fwrite(mfea(i).radvet.vet, sizeof(double), nrad,   arq);
    fwrite(mfea(i).cirvet.vet, sizeof(double), 2*ncir, arq);
  }
  fclose(arq);
}

void le(MFEA& mfea, string nome, int& tnl)
{ FILE* arq=fopen(nome.c_str(),"rb");

  char st[9]; fread(st,1,8,arq); st[7]=0;
  if (string(st)!="IMGFEA*") erro("Erro le(MFEA) - Formato incorreto");

  int nl,nc;
  int nrad,ncir;
  fread(&nl,sizeof(int),1,arq);
  fread(&nc,sizeof(int),1,arq);
  fread(&nrad,sizeof(int),1,arq);
  fread(&ncir,sizeof(int),1,arq);
  fread(&tnl,sizeof(int),1,arq);

  mfea.resize(nl,nc);
  for (int i=0; i<mfea.n(); i++) { 
    mfea.atf(i).rotvet.resize(nrad-1);
    mfea.atf(i).radvet.resize(nrad);
    mfea.atf(i).cirvet.resize(2*ncir);
  }

  for (int i=0; i<mfea.n(); i++) {
    fread(&(mfea(i).rotmap),  sizeof(double), 1,      arq);
    fread(mfea(i).rotvet.vet, sizeof(double), nrad-1, arq);
    fread(mfea(i).radvet.vet, sizeof(double), nrad,   arq);
    fread(mfea(i).cirvet.vet, sizeof(double), 2*ncir, arq);
  }
  fclose(arq);
}

void features(int argc, char** argv)
{ if (argc!=6) {
    printf("Features: Compute rotation-discriminating and -invariant features.\n");
    printf("Features analyze.tga tnl saida.img nrad ncir\n");
    printf("  nrad>=1; nl deve ser impar\n");
    erro("Error: Invalid number of arguments");
  }

  IMGGRY an; pvLe(an,argv[1]);

  int tnl;
  if (sscanf(argv[2],"%d",&tnl)!=1) erro("Erro: Leitura tnl");
  if (tnl<2 || tnl%2!=1) erro("Erro: tnl deve ser impar positivo");
  int raio=tnl/2;

  //IMGGRY anv(an.nl()-tnl+1,an.nc()-tnl+1);
  //for (int l=0; l<anv.nl(); l++)
  //  for (int c=0; c<anv.nc(); c++)
  //    anv(l,c)=an(l+raio,c+raio);

  int nrad;
  if (sscanf(argv[4],"%d",&nrad)!=1) erro("Erro: Leitura nrad");
  if (nrad<1) erro("Erro: nrad<1");

  int ncir;
  if (sscanf(argv[5],"%d",&ncir)!=1) erro("Erro: Leitura ncir");
  if (ncir<0) erro("Erro: ncir<0");

  PESOTOTAL pesototal(nrad,ncir);

  I3DCPX qux; // nrad+ncir, tnl, tnl
  geraqux(qux,nrad,ncir,tnl);

  I3DCPX anc;
  convqux(an,qux,anc);

  IMGGRY valido;
  MFEA anf; // 2*nrad+2*ncir, anv.nl(), anv.nc()
  gerafea(anc,anf,valido,nrad,ncir,tnl);
  //printf("%d %d\n",anf.rotmap.nl(),anf.rotmap.nc());

  imp(anf,argv[3],tnl);
}

void feadist(int argc, char** argv)
{ if (argc!=5) {
    printf("FeaDist: Matching AxT com features pre-calculados. Saida IMGFLT.\n");
    printf("FeaDist analyze.tga features.img template.tga saida.img\n");
    erro("Error: Invalid number of arguments");
  }

  IMGGRY an; pvLe(an,argv[1]);
  int tnl;
  MFEA anf; le(anf,argv[2],tnl);
  int nrad=anf(0).radvet.n();
  int ncir=anf(0).cirvet.n()/2;

  IMGGRY te; pvLe(te,argv[3]);
  if (te.nl()!=te.nc() || te.nl()%2!=1) erro("Erro: template deve ser quadrado de lados impares");
  if (tnl!=te.nl()) erro("Erro: lado do template diferente do tnl do features.img");
  int raio=tnl/2;

  //IMGGRY anv(an.nl()-tnl+1,an.nc()-tnl+1);
  //for (int l=0; l<anv.nl(); l++)
  //  for (int c=0; c<anv.nc(); c++)
  //    anv(l,c)=an(l+raio,c+raio);

  PESOTOTAL pesototal(nrad,ncir);

  I3DCPX qux; // nrad+ncir, tnl, tnl
  geraqux(qux,nrad,ncir,tnl);

  I3DCPX tec;
  convqux(te,qux,tec);

  IMGGRY valido;

  MFEA mtef; // 2*nrad+2*ncir, quv.nl(), quv.nc()
  gerafea(tec,mtef,valido,nrad,ncir,tnl);
  FEA tef=mtef(0); // feature no centro do template
  //printf("%d %d\n",tef.rotmap.nl(),tef.rotmap.nc());

  IMGFLT mat;
  procura(anf,tef,mat,pesototal);
  mat=aumcanvas(mat,an.nl(),an.nc(),1.0);
  pvImp(mat,argv[4]);

  //mat=normaliza(mat); // Tem que mudar normaliza
  //IMGGRY s(an.nl(),an.nc(),255);
  //for (int l=0; l<mat.nl(); l++)
  //  for (int c=0; c<mat.nc(); c++)
  //    s(l+raio,c+raio)=F2G(mat(l,c));
  //pvImp(s,argv[4]);
}

void matdist(int argc, char** argv)
{ if (argc!=6) {
    printf("MatDist: (TeMatch) Matching AxT. Saida grayscale matching distance.\n");
    printf("MatDist analyze.tga template.tga saida.tga nrad ncir\n");
    printf("  nrad>=1\n");
    erro("Error: Invalid number of arguments");
  }

  IMGGRY an; pvLe(an,argv[1]);

  IMGGRY te; pvLe(te,argv[2]);
  if (te.nl()!=te.nc() || te.nl()%2!=1) erro("Erro: template deve ser quadrado de lados impares");
  int tnl=te.nl();
  int raio=tnl/2;

  //IMGGRY anv(an.nl()-tnl+1,an.nc()-tnl+1);
  //for (int l=0; l<anv.nl(); l++)
  //  for (int c=0; c<anv.nc(); c++)
  //    anv(l,c)=an(l+raio,c+raio);

  int nrad;
  if (sscanf(argv[4],"%d",&nrad)!=1) erro("Erro: Leitura nrad");
  if (nrad<1) erro("Erro: nrad<1");

  int ncir;
  if (sscanf(argv[5],"%d",&ncir)!=1) erro("Erro: Leitura ncir");
  if (ncir<0) erro("Erro: ncir<0");

  PESOTOTAL pesototal(nrad,ncir);

  int t1=centseg();

  I3DCPX qux; // nrad+ncir, tnl, tnl
  geraqux(qux,nrad,ncir,tnl);

  I3DCPX tec;
  convqux(te,qux,tec);

  I3DCPX anc;
  convqux(an,qux,anc);

  IMGGRY valido;

  MFEA mtef; // 2*nrad+2*ncir, quv.nl(), quv.nc()
  gerafea(tec,mtef,valido,nrad,ncir,tnl);
  FEA tef=mtef(0); // feature no centro do template
  //printf("%d %d\n",tef.rotmap.nl(),tef.rotmap.nc());

  MFEA anf; // 2*nrad+2*ncir, quv.nl(), quv.nc()
  gerafea(anc,anf,valido,nrad,ncir,tnl);
  //printf("%d %d\n",anf.rotmap.nl(),anf.rotmap.nc());

  int t2=centseg();

  IMGFLT mat;
  procura(anf,tef,mat,pesototal);
  pvImp(mat,argv[3]);

  //vector<PONTOI2> v=kmin(f,1,1);

  //mat=normaliza(mat); // Tem que mudar
  //IMGGRY s(an.nl(),an.nc(),255);
  //for (int l=0; l<mat.nl(); l++)
  //  for (int c=0; c<mat.nc(); c++)
  //    s(l+raio,c+raio)=F2G(mat(l,c));
  //pvImp(s,argv[3]);

  int t3=centseg();
  printf("Pre-processamento=%d processamento=%d\n",t2-t1,t3-t2);
}

double ncc(IMGGRY& an, IMGGRY& te, int lo, int co, double angulo)
// normalize cross correlation 
{ te.backg()=128;
  int raio=te.nl()/2;

  IMGGRY rot=rotacao(te,rad2deg(angulo));
  int conta=0;
  for (int x=rot.minx(); x<=rot.maxx(); x++)
    for (int y=rot.miny(); y<=rot.maxy(); y++) {
      int r=teto(sqrt(double(x*x+y*y)));
      if (r<=raio) {
        if (rot.atC(x,y)==0) rot.atC(x,y)=1;
        conta++;
      } else {
        rot.atC(x,y)=0;
      }
    }
  //pvMostra(rot);

  VETOR<double> Y(conta);
  VETOR<double> X(conta);

  int i=0;
  for (int l=0; l<rot.nl(); l++)
    for (int c=0; c<rot.nc(); c++) {
      if (rot(l,c)>0) {
        Y(i)=an(lo+l,co+c);
        X(i)=rot(l,c);
        i++;
      }
    }
  double coef,beta,gama;
  coefcorr(Y,X,0.001,coef,beta,gama);
  return coef;
}

double mncc(IMGGRY& an, IMGGRY& te, int lo, int co, double& angrad)
// Esta rotina precisa ser otimizada
{ te.backg()=128;
  int raio=te.nl()/2;
  double maxcoef=-infinito; double maxang=angrad;
  double passo=deg2rad(9.0); int npasso=1; int maxind=0;
  //printf("angrad inic de mncc=%g\n",angrad);

  for (int i=-npasso; i<=npasso; i++) {
    IMGGRY rot=rotacao(te,rad2deg(angrad+i*passo));
    int conta=0;
    for (int x=rot.minx(); x<=rot.maxx(); x++)
      for (int y=rot.miny(); y<=rot.maxy(); y++) {
        int r=teto(sqrt(double(x*x+y*y)));
        if (r<=raio && modulo(x,3)==0 && modulo(y,3)==0) {
          if (rot.atC(x,y)==0) rot.atC(x,y)=1;
          conta++;
        } else {
          rot.atC(x,y)=0;
        }
      }
    //pvMostra(rot);

    VETOR<double> Y(conta);
    VETOR<double> X(conta);

    int j=0;
    for (int x=rot.minx(); x<=rot.maxx(); x++)
      for (int y=rot.miny(); y<=rot.maxy(); y++) {
        if (rot.atC(x,y)>0) {
          double r=sqrt(x*x+y*y);
          Y(j)=sqrt(raio-r+1)*an(lo-y+raio,co+x+raio);
          X(j)=sqrt(raio-r+1)*rot.atC(x,y);
          j++;
        }
      }

    double coef,beta,gama;
    coefcorr(Y,X,0.001,coef,beta,gama);
    if (coef>maxcoef) { maxcoef=coef; maxang=angrad+i*passo; maxind=i; }
  }
  angrad=maxang;
  //printf("angrad fim de mncc=%g passo=%g maxind=%d\n",angrad,passo,maxind);
  return maxcoef;
}

void ordena(VETOR<double>& a, vector<PONTOI2>& v, int l, int r)
{ if (!(l<=r)) return;
  int i=l; int j=r;
  double x=a((l+r)/2);
  do {
    while (a(i)<x) i++;
    while (x<a(j)) j--;
    if (i<=j) { swap(a(i),a(j)); swap(v[i],v[j]); i++; j--; }
  } while (i<=j);
  if (l<j) ordena(a,v,l,j);
  if (i<r) ordena(a,v,i,r);
}

void mordena(VETOR<double>& a, vector<PONTOI2>& v, VETOR<double>& vang, int l, int r)
{ if (!(l<=r)) return;
  int i=l; int j=r;
  double x=a((l+r)/2);
  do {
    while (a(i)<x) i++;
    while (x<a(j)) j--;
    if (i<=j) { swap(a(i),a(j)); swap(v[i],v[j]); swap(vang(i),vang(j)); i++; j--; }
  } while (i<=j);
  if (l<j) mordena(a,v,vang,l,j);
  if (i<r) mordena(a,v,vang,i,r);
}

void mordena2(VETOR<double>& a, vector<PONTOI2>& v, VETOR<double>& vang, VETOR<int>& vesc, int l, int r)
{ if (!(l<=r)) return;
  int i=l; int j=r;
  double x=a((l+r)/2);
  do {
    while (a(i)<x) i++;
    while (x<a(j)) j--;
    if (i<=j) { swap(a(i),a(j)); swap(v[i],v[j]); swap(vang(i),vang(j)); swap(vesc(i),vesc(j)); i++; j--; }
  } while (i<=j);
  if (l<j) mordena2(a,v,vang,vesc,l,j);
  if (i<r) mordena2(a,v,vang,vesc,i,r);
}

void mordena3(VETOR<double>& a, vector<PONTOI2>& v, VETOR<double>& vang, VETOR<int>& vesc, VETOR<int>& vq, int l, int r)
{ if (!(l<=r)) return;
  int i=l; int j=r;
  double x=a((l+r)/2);
  do {
    while (a(i)<x) i++;
    while (x<a(j)) j--;
    if (i<=j) { swap(a(i),a(j)); swap(v[i],v[j]); swap(vang(i),vang(j)); swap(vesc(i),vesc(j)); swap(vq(i),vq(j)); i++; j--; }
  } while (i<=j);
  if (l<j) mordena3(a,v,vang,vesc,vq,l,j);
  if (i<r) mordena3(a,v,vang,vesc,vq,i,r);
}

void gerate(int argc, char** argv)
{ if (argc!=10) {
    printf("GeraTe: Gera templates a partir de query\n");
    printf("GeraTe query.tga template nrad ncir tnl precisao nte separacao v/s\n");
    printf("  template sem extensao\n");
    printf("  nrad>=1\n");
    printf("  precisao entre 0 e 1\n");
    printf("  Se nte==1, deve especificar template.tga (com extensao)\n");
    erro("Error: Invalid number of arguments");
  }

  IMGGRY qu; pvLe(qu,argv[1]);

  int nrad;
  if (sscanf(argv[3],"%d",&nrad)!=1) erro("Erro: Leitura nrad");
  if (nrad<1) erro("Erro: nrad<1");

  int ncir;
  if (sscanf(argv[4],"%d",&ncir)!=1) erro("Erro: Leitura ncir");
  if (ncir<0) erro("Erro: ncir<0");

  PESOTOTAL pesototal(nrad,ncir);

  int tnl;
  if (sscanf(argv[5],"%d",&tnl)!=1) erro("Erro: Leitura tnl");
  if (tnl%2==0 || tnl<=0) erro("Erro: tnl deve ser impar");
  int raio=tnl/2;

  IMGGRY quv(qu.nl()-tnl+1,qu.nc()-tnl+1);
  for (int l=0; l<quv.nl(); l++)
    for (int c=0; c<quv.nc(); c++)
      quv(l,c)=qu(l+raio,c+raio);

  double precisao;
  if (sscanf(argv[6],"%lf",&precisao)!=1) erro("Erro: Leitura precisao");
  if (precisao<0.0 || 1.0<precisao) erro("Erro: Precisao fora de [0,1]");

  int nte;
  if (sscanf(argv[7],"%d",&nte)!=1) erro("Erro: Leitura nte");
  if (nte<1) erro("Erro: nte<1");

  int separacao;
  if (sscanf(argv[8],"%d",&separacao)!=1) erro("Erro: Leitura separacao");
  if (separacao<1) erro("Erro: separacao<1");

  char verbose=tolower(argv[9][0]);
  if (verbose!='v' && verbose!='s') erro("Erro: Verbose deve ser v/s");

  I3DCPX qux; // nrad+ncir, tnl, tnl
  geraqux(qux,nrad,ncir,tnl);

  I3DCPX quc;
  convqux(qu,qux,quc);

  MFEA quf; // 2*nrad+2*ncir, quv.nl(), quv.nc()
  IMGGRY valido;
  gerafea(quc,quf,valido,nrad,ncir,tnl);
  //pvMostra(valido,"valido minampl");

  stablecand(quf,valido,precisao,pesototal);
  //pvMostra(valido,"valido stable");
  //pvImp(valido,"stable.tga");

  //if (nerosao>1) valido = dilat(valido,nerosao,nerosao);
  //valido = findce(valido);
  //pvMostra(valido,"valido apos findce",1);

  IMGFLT f=distgeo(valido);
  //pvMostra(IMGGRY(normaliza(f)),"valido distgeo");

  vector<PONTOI2> v=kmax(f,nte,separacao);
  //for (unsigned i=0; i<v.size(); i++) cout << v[i];
  //cout << endl;
  if (v.size()==0) erro("Erro: Sem ponto estavel");

  if (verbose=='v') {
    IMGCOR b=qu;
    for (int l=0; l<quv.nl(); l++)
      for (int c=0; c<quv.nc(); c++)
        if (valido(l,c)==0) b(l+raio,c+raio)=COR(0,0,b(l+raio,c+raio).b());
    for (unsigned i=0; i<v.size(); i++) b(v[i].l()+raio,v[i].c()+raio)=COR(255,0,0);
    pvImp(b,"valido.tga");
  }

  if (nte>1) {
    IMGGRY sai(tnl,tnl);
    for (unsigned s=0; s<v.size(); s++) {
      int l=v[s].l(); int c=v[s].c();
      for (int l2=-raio; l2<=raio; l2++)
        for (int c2=-raio; c2<=raio; c2++)
          sai.atc(l2,c2)=qu(l+l2+raio,c+c2+raio);
      char st[256];
      sprintf(st,"%s%02d.tga",argv[2],s);
      pvImp(sai,st);
    }
  } else {
    IMGGRY sai(tnl,tnl);
    int l=v[0].l(); int c=v[0].c();
    for (int l2=-raio; l2<=raio; l2++)
      for (int c2=-raio; c2<=raio; c2++)
        sai.atc(l2,c2)=qu(l+l2+raio,c+c2+raio);
    pvImp(sai,argv[2]);
  }
}

//<<<<<<<<<<<<<<<<<<< teruim <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
PONTOI2 candruim(MFEA& mfea, PESOTOTAL& pesototal, char modo)
{ double maxruim=-infinito;
  PONTOI2 p(0,0);
  for (int l=1; l<mfea.nl()-1; l++)
    for (int c=1; c<mfea.nc()-1; c++) {

      double maxdif1=-infinito;
      for (int l2=-1; l2<=1; l2++)
        for (int c2=-1; c2<=1; c2++) {
          double t=difangrad(mfea.atf(l,c).rotmap,mfea.atf(l+l2,c+c2).rotmap)/M_PI;
          if (t>maxdif1) maxdif1=t;
        }

      double maxdif2=-infinito;
      for (int l2=-1; l2<=1; l2++)
        for (int c2=-1; c2<=1; c2++) {
          double t=diffea(mfea.atf(l,c),mfea.atf(l+l2,c+c2),pesototal);
          if (t>maxdif2) maxdif2=t;
        }

      double maxdif=0.0;
      if (modo=='o') {
        maxdif=maxdif1;
      } else if (modo=='f') {
        maxdif=maxdif2;
      } else if (modo=='a') {
        maxdif=(maxdif1+maxdif2)/2.0;
      } else if (modo=='m') {
        maxdif=min(maxdif1,maxdif2);
      }
      if (maxdif>maxruim) { maxruim=maxdif; p=PONTOI2(l,c); }
    }
  printf("maxruim=%f\n",maxruim);
  return p;
}

void teruim(int argc, char** argv)
{ if (argc!=7) {
    printf("TeRuim: Gera propositalmente um template ruim a partir de query\n");
    printf("TeRuim query.tga template.tga nrad ncir tnl o/f/a/m\n");
    printf("  nrad>=1\n");
    printf("  modo=o orientacao; modo=f feature; modo=a average; modo=m minimo\n");
    erro("Error: Invalid number of arguments");
  }

  IMGGRY qu; pvLe(qu,argv[1]);

  int nrad;
  if (sscanf(argv[3],"%d",&nrad)!=1) erro("Erro: Leitura nrad");
  if (nrad<1) erro("Erro: nrad<1");

  int ncir;
  if (sscanf(argv[4],"%d",&ncir)!=1) erro("Erro: Leitura ncir");
  if (ncir<0) erro("Erro: ncir<0");

  PESOTOTAL pesototal(nrad,ncir);

  int tnl;
  if (sscanf(argv[5],"%d",&tnl)!=1) erro("Erro: Leitura tnl");
  if (tnl%2==0 || tnl<=0) erro("Erro: tnl deve ser impar");
  int raio=tnl/2;

  char modo=tolower(argv[6][0]);
  if (modo!='o' && modo!='f' && modo!='a' && modo!='m') erro("Erro: Modo invalido");

  IMGGRY quv(qu.nl()-tnl+1,qu.nc()-tnl+1);
  for (int l=0; l<quv.nl(); l++)
    for (int c=0; c<quv.nc(); c++)
      quv(l,c)=qu(l+raio,c+raio);

  I3DCPX qux; // nrad+ncir, tnl, tnl
  geraqux(qux,nrad,ncir,tnl);

  I3DCPX quc;
  convqux(qu,qux,quc);

  MFEA quf; // 2*nrad+2*ncir, quv.nl(), quv.nc()
  IMGGRY valido;
  gerafea(quc,quf,valido,nrad,ncir,tnl);
  //pvMostra(valido,"valido minampl");

  PONTOI2 p=candruim(quf,pesototal,modo);

  IMGGRY sai(tnl,tnl);
  int l=p.l(); int c=p.c();
  for (int l2=-raio; l2<=raio; l2++)
    for (int c2=-raio; c2<=raio; c2++)
      sai.atc(l2,c2)=qu(l+l2+raio,c+c2+raio);
  pvImp(sai,argv[2]);
}

void mtefi1(int argc, char** argv)
{
  if (argc<9) {
    printf("MTefi1: Matching de A e varios Ts com verificacao NCC. A contem 1 unico T.\n");
    printf("MTefi1 A.tga P.tga nrad ncir prec detec ncand T1.tga T2.tga ...\n");
    printf("  nrad>=1; ncand<=1 => nao faz NCC\n");
    printf("  Ex: a.tga p.tga 3 3 0.15 1.2 40 t*.tga\n");
    printf("  Todos Ts devem ter mesmo nl=nc impar\n");
    erro("Error: Invalid number of arguments");
  }

  IMGGRY an; pvLe(an,argv[1]);

  IMGGRY te; pvLe(te,argv[8]);
  if (te.nl()!=te.nc() || te.nl()%2!=1) erro("Erro: template deve ser quadrado de lados impares");
  int tnl=te.nl();
  int raio=tnl/2;

  int t1=centseg();

  int nrad;
  if (sscanf(argv[3],"%d",&nrad)!=1) erro("Erro: Leitura nrad");
  if (nrad<1) erro("Erro: nrad<1");

  int ncir;
  if (sscanf(argv[4],"%d",&ncir)!=1) erro("Erro: Leitura ncir");
  if (ncir<0) erro("Erro: ncir<0");

  double precisao;
  if (sscanf(argv[5],"%lf",&precisao)!=1) erro("Erro: Leitura precisao");
  if (precisao<0.0 || 1.0<precisao) erro("Erro: Precisao fora de [0,1]");

  double detec;
  if (sscanf(argv[6],"%lf",&detec)!=1) erro("Erro: Leitura detec");
  if (detec<1.0) erro("Erro: detec<1.0");

  int ncand;
  if (sscanf(argv[7],"%d",&ncand)!=1) erro("Erro: Leitura ncand");
  if (ncand<0) erro("Erro: ncand<0");

  IMGGRY anv(an.nl()-tnl+1,an.nc()-tnl+1);
  for (int l=0; l<anv.nl(); l++)
    for (int c=0; c<anv.nc(); c++) {
      anv(l,c)=an(l+raio,c+raio);
  }

  double deteccao=detec*precisao;
  PESOTOTAL pesototal(nrad,ncir);

  I3DCPX qux; // nrad+ncir, tnl, tnl
  geraqux(qux,nrad,ncir,tnl);

  I3DCPX anc;
  convqux(an,qux,anc);

  int t2=centseg();

  MFEA anf; // 2*nrad+2*ncir, quv.nl(), quv.nc()
  IMGGRY valido;
  gerafea(anc,anf,valido,nrad,ncir,tnl);
  //printf("%d %d\n",anf.rotmap.nl(),anf.rotmap.nc());

  IMGCOR sai=anv;
  int t3=centseg();

  for (int i=8; i<argc; i++) {
    pvLe(te,argv[i]);
    if (te.nl()!=te.nc() || te.nl()%2!=1 || te.nl()!=tnl) erro("Erro: template com lado diferente");

    I3DCPX tec;
    convqux(te,qux,tec);

    MFEA mtef; // 2*nrad+2*ncir, quv.nl(), quv.nc()
    gerafea(tec,mtef,valido,nrad,ncir,tnl);
    FEA tef=mtef(0);
    //printf("%d %d\n",tef.rotmap.nl(),tef.rotmap.nc());

    IMGFLT mat;
    // procura1(anf,tef,mat,deteccao,pesototal);
    // Acho que nao pode usar procura1 aqui, pois vai procurar ncand matchings
    procura(anf,tef,mat,pesototal); 


    PONTOI2 p;
    double angrad;
    if (ncand>1) {
      vector<PONTOI2> v=kargmin(mat,ncand);
      if (int(v.size())!=ncand)
        printf("Numero de candidatos procurados=%d achados=%d\n",ncand,v.size());
      if (v.size()==0) erro("Erro: Nao achei candidato");

      //for (unsigned i=0; i<v.size(); i++)
      //  cruz(sai,v[i].l(),v[i].c(),3,COR(0,0,128));
      //pvImp(sai,argv[3]);

      VETOR<double> match(v.size());
      VETOR<double> vang(v.size());
      for (unsigned i=0; i<v.size(); i++) {
        int l=v[i].l(); int c=v[i].c();
        double angrad2=anf.atf(l,c).rotmap-tef.rotmap; // em radianos
        match(i)=mncc(an,te,l,c,angrad2);
        vang(i)=angrad2;
      }
      //cout << "match: " << match;
      mordena(match,v,vang,0,match.n()-1);
      //cout << "match: " << match;
      //for (unsigned i=0; i<v.size(); i++) cout << v[i] << "  ";
      p=v[v.size()-1];
      angrad=-vang(v.size()-1);
    } else {
      p=argmin(mat);
      int l=p.l(); int c=p.c();
      angrad=-anf.atf(l,c).rotmap+tef.rotmap; // em radianos
    }
    relogioretangular2(sai,p.l(),p.c(),angrad,tnl,tnl,i-7);
  }
  int t4=centseg();

  pvImp(sai,argv[2]);
  printf("Conv_a=%d gerafea_a=%d processa_t=%d total=%d\n",t2-t1,t3-t2,t4-t3,t4-t1);
}

void mvhfi1(int argc, char** argv)
{
  if (argc<12) {
    printf("MVHfi1: Matching A x varios Qs com VHough. Gera Ts. Cd Q aparece 1 vez em A\n");
    printf("MVHfi1 A.tga P.tga nrad ncir tnl prec nte detec ncand separ Q*.tga\n");
    printf("Ex: nrad=3 ncir=3 tnl=21 prec=0.1 nte=8 detec=1.2 ncand=10 separ=11 Q*.tga\n");
    printf("  Prec e' usado na geracao de candidatos estaveis\n");
    printf("  Deteccao=prec*detec usado na procura1\n");
    printf("  Para cada T, ncand pixels com menores distancias entram no Hough\n");
    erro("Error: Invalid number of arguments");
  }

  int t1=centseg();
  IMGGRY an; pvLe(an,argv[1]);

  int nrad;
  if (sscanf(argv[3],"%d",&nrad)!=1) erro("Erro: Leitura nrad");
  if (nrad<1) erro("Erro: nrad<1");

  int ncir=3;
  if (sscanf(argv[4],"%d",&ncir)!=1) erro("Erro: Leitura ncir");
  if (ncir<0) erro("Erro: ncir<0");

  int tnl;
  if (sscanf(argv[5],"%d",&tnl)!=1) erro("Erro: Leitura tnl");
  if (tnl%2==0 || tnl<=0) erro("Erro: tnl deve ser impar");

  double precisao;
  if (sscanf(argv[6],"%lf",&precisao)!=1) erro("Erro: Leitura precisao");
  if (precisao<0.0 || 1.0<precisao) erro("Erro: Precisao fora de [0,1]");

  int nte;
  if (sscanf(argv[7],"%d",&nte)!=1) erro("Erro: Leitura nte");
  if (nte<1) erro("Erro: nte<1");

  double detec;
  if (sscanf(argv[8],"%lf",&detec)!=1) erro("Erro: Leitura detec");
  if (detec<=1.0) erro("Erro: detec<=1.0");
  double deteccao=detec*precisao;

  int ncand;
  if (sscanf(argv[9],"%d",&ncand)!=1) erro("Erro: Leitura ncand");
  if (ncand<1) erro("Erro: ncand<1");

  int separacao;
  if (sscanf(argv[10],"%d",&separacao)!=1) erro("Erro: Leitura separacao");
  if (separacao<0 || tnl<separacao) erro("Erro: separacao<0 || tnl<separacao");

  int raio=tnl/2;

  IMGGRY anv(an.nl()-tnl+1,an.nc()-tnl+1);
  for (int l=0; l<anv.nl(); l++)
    for (int c=0; c<anv.nc(); c++) {
      anv(l,c)=an(l+raio,c+raio);
  }

  PESOTOTAL pesototal(nrad,ncir);

  I3DCPX qux; // nrad+ncir, tnl, tnl
  geraqux(qux,nrad,ncir,tnl);

  I3DCPX anc;
  convqux(an,qux,anc);

  //IMGGRY valido;
  MFEA anf; // 2*nrad+2*ncir, quv.nl(), quv.nc()
  IMGGRY valido;
  gerafea(anc,anf,valido,nrad,ncir,tnl);

  IMGGRY qu;
  IMGCOR sai=an;
  int arginic=11;

  int t2=centseg();

  int tprocura=0; int though=0;
  for (int q=arginic; q<argc; q++) {
    int tl3=centseg();
    pvLe(qu,argv[q]);
    printf("Processando %s\n",argv[q]);
    //pvMostra(qu);

    IMGGRY quv(qu.nl()-tnl+1,qu.nc()-tnl+1);
    for (int l=0; l<quv.nl(); l++)
      for (int c=0; c<quv.nc(); c++)
        quv(l,c)=qu(l+raio,c+raio);

    I3DCPX quc;
    convqux(qu,qux,quc);

    MFEA quf; // 2*nrad+2*ncir, quv.nl(), quv.nc()
    gerafea(quc,quf,valido,nrad,ncir,tnl);

    stablecand(quf,valido,precisao,pesototal);

    IMGFLT f=distgeo(valido);
    //pvMostra(IMGGRY(normaliza(f)),"valido distgeo",1);
    vector<PONTOI2> vt=kmax(f,nte,separacao);
    if (nte!=int(vt.size())) printf("Warning: query argv[%d]. #templ procurados=%d achados=%d\n",q,nte,vt.size());
    if (vt.size()==0) erro("Erro: Sem ponto estavel");

    IMGFLT mat;
    FEA tef;

    //printf("ncand=%d\n",ncand);

    MATRIZ<PONTOI2> pos(nte, ncand,PONTOI2(maxint,maxint));
    IMGCPX  ana(nte, ncand,CPX(0,0));
    IMGFLT  hou(nte, ncand,-infinito);

    int tl4=centseg(); tprocura += (tl4-tl3);

    double passo=max(tnl/5,3); // tolerancia
    for (unsigned i=0; i<vt.size(); i++) {
      tef=quf(vt[i].l(),vt[i].c());
      int dl=quv.lc()-vt[i].l(); int dc=quv.cc()-vt[i].c();
      //procura1(anf,tef,mat,deteccao,pesototal);
      // Acho que nao pode usar procura1 aqui, pois vai procurar ncand matchings
      procura(anf,tef,mat,pesototal); 

      vector<PONTOI2> vh=kargmin(mat,ncand);
      if (ncand!=int(vh.size())) 
         printf("Warning: query argv[%d] templ=%d: #cand procurados=%d achados=%d\n",q,i,ncand,vh.size());

      for (unsigned j=0; j<vh.size(); j++) {
        int l=vh[j].l(); int c=vh[j].c();
        double peso=1.0-mat.atf(l,c);
        double angulo=-anf.atf(l,c).rotmap+tef.rotmap;
        double co=cos(angulo); double si=sin(angulo);
        int ndl=arredonda( co*dl+si*dc);
        int ndc=arredonda(-si*dl+co*dc);
        hou.atf(i,j)=peso;
        ana.atf(i,j)=CPX(peso*co,peso*si);
        pos.atf(i,j)=PONTOI2(l+ndl,c+ndc);
      }
    }

    double maxsomapeso=-infinito;
    PONTO2 maxmedia(0.0,0.0);
    CPX maxana(0.0,0.0);
    int maxconta=0;
    for (int c=0; c<pos.nc(); c++) {
      for (int l=0; l<pos.nl(); l++) { // Tem que fazer loop c/l para pegar melhor matching
        double somapeso=hou.atf(l,c);
        PONTO2 somaponto=somapeso*PONTO2(pos.atf(l,c));
        CPX somaana=ana.atf(l,c);
        int conta=1;
        for (int l2=0; l2<pos.nl(); l2++) {
          if (l==l2) break;
          for (int c2=0; c2<pos.nc(); c2++) {
            if (distancia(pos.atf(l,c),pos.atf(l2,c2))<passo) {
              somapeso += hou.atf(l2,c2);
              somaponto = somaponto + hou.atf(l2,c2)*PONTO2(pos.atf(l2,c2));
              somaana += ana.atf(l2,c2); // Nao esta multiplicando peso
              conta++;
              break;
            }
          }
        }
        if (somapeso>maxsomapeso) {
          maxsomapeso=somapeso;
          maxmedia=(1.0/somapeso)*somaponto;
          maxana=somaana;
          maxconta=conta;
          if (maxconta==pos.nl()) goto saida;
        }
      }
    }
    saida:

    double angulo=arg(maxana);
    int l=arredonda(maxmedia.l())+raio; int c=arredonda(maxmedia.c())+raio;
    relogioretangular2(sai,l,c,angulo,qu.nl(),qu.nc(),q-arginic+1);
    int tl5=centseg(); though += (tl5-tl4);
  }
  pvImp(sai,argv[2]);
  int t3=centseg();
  printf("Conv_fea_a=%d proc_hough=%d total=%d procura=%d hough=%d\n",t2-t1,t3-t2,t3-t1,tprocura,though);
}

void mstefi1(int argc, char** argv)
{
  if (argc<13) {
    printf("MSTefi1: Scale invariant MTefi com teste NCC. Cd Q aparece 1 vez em A\n");
    printf("MSTefi1 A.tga P.tga nrad ncir prec detec ncand escinic escfim nesc tnl Q*.tga\n");
    printf("  nrad>=1; ncand<=1 => nao faz NCC\n");
    printf("  Ex: a.tga p.tga 4 4 0.1 1.6 10 0.7 1.4 8 31 q*.tga\n");
    erro("Error: Invalid number of arguments");
  }

  int t1=centseg();
  IMGGRY an; pvLe(an,argv[1]);

  int nrad;
  if (sscanf(argv[3],"%d",&nrad)!=1) erro("Erro: Leitura nrad");
  if (nrad<1) erro("Erro: nrad<1");

  int ncir;
  if (sscanf(argv[4],"%d",&ncir)!=1) erro("Erro: Leitura ncir");
  if (ncir<0) erro("Erro: ncir<0");

  double precisao;
  if (sscanf(argv[5],"%lf",&precisao)!=1) erro("Erro: Leitura precisao");
  if (precisao<0.0 || 1.0<precisao) erro("Erro: Precisao fora de [0,1]");

  double detec;
  if (sscanf(argv[6],"%lf",&detec)!=1) erro("Erro: Leitura detec");
  if (detec<1.0) erro("Erro: detec<1.0");

  int ncand;
  if (sscanf(argv[7],"%d",&ncand)!=1) erro("Erro: Leitura ncand");
  if (ncand<0) erro("Erro: ncand<0");

  double escinic;
  if (sscanf(argv[8],"%lf",&escinic)!=1) erro("Erro: Leitura escinic");
  if (escinic<0.0) erro("Erro: escinic<0.0");

  double escfim;
  if (sscanf(argv[9],"%lf",&escfim)!=1) erro("Erro: Leitura escfim");
  if (escfim<escinic) erro("Erro: escfim<escinic");

  int nesc;
  if (sscanf(argv[10],"%d",&nesc)!=1) erro("Erro: Leitura nesc");
  if (nesc<1) erro("Erro: nesc<1");
  double passoesc=exp(log(escfim/escinic)/(nesc-1));

  int tnl;
  if (sscanf(argv[11],"%d",&tnl)!=1) erro("Erro: Leitura tnl");
  if (tnl%2==0 || tnl<=0) erro("Erro: tnl deve ser impar");
  int raio=tnl/2;

  IMGGRY anv(an.nl()-tnl+1,an.nc()-tnl+1);
  for (int l=0; l<anv.nl(); l++)
    for (int c=0; c<anv.nc(); c++) {
      anv(l,c)=an(l+raio,c+raio);
  }

  double deteccao=detec*precisao;
  PESOTOTAL pesototal(nrad,ncir);

  I3DCPX qux; // nrad+ncir, tnl, tnl
  geraqux(qux,nrad,ncir,tnl);

  I3DCPX anc;
  convqux(an,qux,anc);

  MFEA anf; // 2*nrad+2*ncir, quv.nl(), quv.nc()
  IMGGRY valido;
  gerafea(anc,anf,valido,nrad,ncir,tnl);
  //printf("%d %d\n",anf.rotmap.nl(),anf.rotmap.nc());

  IMGGRY qu;
  IMGCOR sai=anv;
  int arginic=12;

  int t2=centseg();
  int tgeraproc=0; int tncc=0; int tprocura=0; int tgerate=0;

  for (int q=arginic; q<argc; q++) {
    int tl1=centseg();

    pvLe(qu,argv[q]);

    // Geracao dos Qs e templates
    VETOR<IMGGRY> qesc(nesc);
    for (int s=0; s<nesc; s++) {
      double esc=escinic*pow(passoesc,s);
      int nl=arredonda(esc*qu.nl());
      int nc=arredonda(esc*qu.nc());
      qesc(s)=zoomnp(qu,nl,nc);
    }

    I3DGRY te(nesc,tnl,tnl);
    IMGFLT zbuf(anv.nl(),anv.nc(),infinito);
    IMGSHT zind(anv.nl(),anv.nc(),-1);
    VETOR<FEA> tef(nesc);
    for (int s=0; s<nesc; s++) {
      //printf("Escala s=%d\n",s);
      int tm1=centseg();

      I3DCPX quc;
      convqux(qesc(s),qux,quc);

      MFEA quf; // 2*nrad+2*ncir, quv.nl(), quv.nc()
      IMGGRY valido;
      gerafea(quc,quf,valido,nrad,ncir,tnl);
      stablecand(quf,valido,precisao,pesototal);
      IMGFLT f=distgeo(valido);
      //pvMostra(IMGGRY(normaliza(f)),"valido distgeo",1);
      PONTOI2 p=argmax(f);
      int l1=p.l(); int c1=p.c();
      //p contem indices em quv e quf
      if (f(l1,c1)==0.0) erro("Erro: Sem ponto estavel");
      for (int l2=0; l2<tnl; l2++)
        for (int c2=0; c2<tnl; c2++)
          te(s,l2,c2)=qesc.atf(s).atf(l1+l2,c1+c2);

      tef(s)=quf(l1,c1);

      int tm2=centseg();
      tgerate += (tm2-tm1);

      IMGFLT mat;
      //procura1(anf,tef(s),mat,deteccao,pesototal);
      //if (mat.n()!=anv.n()) erro("Erro");
      // Acho que nao pode usar procura1 aqui, pois vai procurar ncand matchings
      procura(anf,tef(s),mat,pesototal); 

      for (int i=0; i<mat.n(); i++)
        if (mat(i)<zbuf(i)) {
          zbuf(i)=mat(i);
          zind(i)=s;
        }
      int tm3=centseg();
      tprocura += (tm3-tm2);
    }
    int tl2=centseg();
    tgeraproc += (tl2-tl1);

    PONTOI2 p;
    double angrad;
    int s;
    if (ncand>1) {
      vector<PONTOI2> v=kargmin(zbuf,ncand);
      //for (unsigned z=0; z<v.size(); z++) cout << v[z];
      //cout << endl;

      //vector<PONTOI2> v=kmin(zbuf,ncand,tnl/2);
      //for (unsigned z=0; z<v.size(); z++) cout << v[z];
      //cout << endl;
      if (ncand!=int(v.size())) printf("Warning: #candidatos procurados=%d achados=%d\n",ncand,v.size());
      if (v.size()==0) erro("Erro: Nao achei candidato");

      VETOR<double> match(v.size());
      VETOR<double> vang(v.size());
      VETOR<int> vesc(v.size());
      for (unsigned i=0; i<v.size(); i++) {
        int l=v[i].l(); int c=v[i].c();
        int s2=vesc(i)=zind(l,c);
        double angrad2=anf.atf(l,c).rotmap-tef(s2).rotmap; // em radianos
        match(i)=mncc(an,te(s2),l,c,angrad2);
        vang(i)=angrad2;
      }
      mordena2(match,v,vang,vesc,0,match.n()-1);
      p=v[v.size()-1];
      angrad=-vang(v.size()-1);
      s=vesc(v.size()-1); // coloquei arredonda por que dava warning.
    } else {
      p=argmin(zbuf);
      int l=p.l(); int c=p.c();
      s=zind(l,c);
      angrad=-anf.atf(l,c).rotmap+tef(s).rotmap; // em radianos
    }
    int tnl2=arredonda(escinic*pow(passoesc,s)*tnl);
    relogioretangular2(sai,p.l(),p.c(),angrad,tnl2,tnl2,q-arginic+1);
    int tl3=centseg();
    tncc += (tl3-tl2);
  }
  pvImp(sai,argv[2]);

  int t3=centseg();
  printf("Conv_fea_a=%d procura=%d total=%d gerate_procura=%d cand_ncc=%d gerate=%d procura=%d\n",t2-t1,t3-t2,t3-t1,tgeraproc,tncc,tgerate,tprocura);
}

void aftefi1(int argc, char** argv)
{ if (argc<13) {
    printf("AfTefi1: Affine invariant MTefi com teste NCC. Cd Q aparece 1 vez em A\n");
    printf("  Deve fornecer Qs distorcidos\n");
    printf("AfTefi1 A.tga P.tga nrad ncir prec detec ncand escinic escfim nesc tnl Q*.tga\n");
    printf("  nrad>=1; ncand<=1 => nao faz NCC\n");
    printf("  Ex: a.tga p.tga 4 4 0.1 1.6 10 0.7 1.4 8 31 q*.tga\n");
    erro("Error: Invalid number of arguments");
  }

  int t1=centseg();
  IMGGRY an; pvLe(an,argv[1]);

  int nrad;
  if (sscanf(argv[3],"%d",&nrad)!=1) erro("Erro: Leitura nrad");
  if (nrad<1) erro("Erro: nrad<1");

  int ncir;
  if (sscanf(argv[4],"%d",&ncir)!=1) erro("Erro: Leitura ncir");
  if (ncir<0) erro("Erro: ncir<0");

  double precisao;
  if (sscanf(argv[5],"%lf",&precisao)!=1) erro("Erro: Leitura precisao");
  if (precisao<0.0 || 1.0<precisao) erro("Erro: Precisao fora de [0,1]");

  double detec;
  if (sscanf(argv[6],"%lf",&detec)!=1) erro("Erro: Leitura detec");
  if (detec<1.0) erro("Erro: detec<1.0");

  int ncand;
  if (sscanf(argv[7],"%d",&ncand)!=1) erro("Erro: Leitura ncand");
  if (ncand<0) erro("Erro: ncand<0");

  double escinic;
  if (sscanf(argv[8],"%lf",&escinic)!=1) erro("Erro: Leitura escinic");
  if (escinic<0.0) erro("Erro: escinic<0.0");

  double escfim;
  if (sscanf(argv[9],"%lf",&escfim)!=1) erro("Erro: Leitura escfim");
  if (escfim<escinic) erro("Erro: escfim<escinic");

  int nesc;
  if (sscanf(argv[10],"%d",&nesc)!=1) erro("Erro: Leitura nesc");
  if (nesc<1) erro("Erro: nesc<1");
  double passoesc=exp(log(escfim/escinic)/(nesc-1));

  int tnl;
  if (sscanf(argv[11],"%d",&tnl)!=1) erro("Erro: Leitura tnl");
  if (tnl%2==0 || tnl<=0) erro("Erro: tnl deve ser impar");
  int raio=tnl/2;

  IMGGRY anv(an.nl()-tnl+1,an.nc()-tnl+1);
  for (int l=0; l<anv.nl(); l++)
    for (int c=0; c<anv.nc(); c++) {
      anv(l,c)=an(l+raio,c+raio);
  }

  double deteccao=detec*precisao;
  PESOTOTAL pesototal(nrad,ncir);

  I3DCPX qux; // nrad+ncir, tnl, tnl
  geraqux(qux,nrad,ncir,tnl);

  I3DCPX anc;
  convqux(an,qux,anc);

  MFEA anf; // 2*nrad+2*ncir, quv.nl(), quv.nc()
  IMGGRY valido;
  gerafea(anc,anf,valido,nrad,ncir,tnl);
  //printf("%d %d\n",anf.rotmap.nl(),anf.rotmap.nc());

  IMGGRY qu;
  IMGCOR sai=anv;
  int arginic=12;
  MATRIZ<IMGGRY> qesc(nesc,argc-arginic);
  MATRIZ<FEA> tef(nesc,argc-arginic);
  MATRIZ<IMGGRY> te(nesc,argc-arginic); te.fill(IMGGRY(tnl,tnl,0));
  int t2=centseg();

  IMGFLT zbuf(anv.nl(),anv.nc(),infinito);
  IMGSHT zinds(anv.nl(),anv.nc(),-1);
  IMGSHT zindq(anv.nl(),anv.nc(),-1);

  for (int q=0; q<qesc.nc(); q++) {
    printf("Processando %s\n",argv[q+arginic]);
    pvLe(qu,argv[q+arginic]);

    // Geracao dos Qs e templates
    for (int s=0; s<nesc; s++) {
      double esc=escinic*pow(passoesc,s);
      int nl=arredonda(esc*qu.nl());
      int nc=arredonda(esc*qu.nc());
      qesc(s,q)=zoomnp(qu,nl,nc);
    }

    for (int s=0; s<nesc; s++) {
      //printf("Escala s=%d\n",s);

      I3DCPX quc;
      convqux(qesc.atf(s,q),qux,quc);

      MFEA quf; // 2*nrad+2*ncir, quv.nl(), quv.nc()
      IMGGRY valido;
      gerafea(quc,quf,valido,nrad,ncir,tnl);
      if (ignora255) Ignora255(qesc.atf(s,q),valido,tnl);
      stablecand(quf,valido,precisao,pesototal);
      IMGFLT f=distgeo(valido);
      //pvMostra(IMGGRY(normaliza(f)),"valido distgeo",1);
      PONTOI2 p=argmax(f);
      int l1=p.l(); int c1=p.c();
      //p contem indices em quv e quf
      if (f(l1,c1)==0.0) { printf("Warning: Sem ponto estavel q=%d s=%d\n",q,s); continue; }
      for (int l2=0; l2<tnl; l2++)
        for (int c2=0; c2<tnl; c2++)
          te(s,q)(l2,c2)=qesc.atf(s,q).atf(l1+l2,c1+c2);

      char st[256];
      sprintf(st,"t_q%02d_s%02d.png",q,s);
      pvImp(te(s,q),st);

      tef(s)=quf(l1,c1);
      IMGFLT mat;
      //procura1(anf,tef(s),mat,deteccao,pesototal);
      //if (mat.n()!=anv.n()) erro("Erro");
      // Acho que nao pode usar procura1 aqui, pois vai procurar ncand matchings
      procura(anf,tef(s),mat,pesototal); 

      for (int i=0; i<mat.n(); i++)
        if (mat(i)<zbuf(i)) {
          zbuf(i)=mat(i);
          zinds(i)=s;
          zindq(i)=q;
        }
    }
  }

  PONTOI2 p;
  double angrad;
  int s=0, q=0;
  if (ncand>1) {
    vector<PONTOI2> v=kargmin(zbuf,ncand);
    if (ncand!=int(v.size())) printf("Warning: #candidatos procurados=%d achados=%d\n",ncand,v.size());
    if (v.size()==0) erro("Erro: Nao achei candidato");

    VETOR<double> match(v.size());
    VETOR<double> vang(v.size());
    VETOR<int> vesc(v.size());
    VETOR<int> vq(v.size());
    for (unsigned i=0; i<v.size(); i++) {
      int l=v[i].l(); int c=v[i].c();
      int s2=vesc(i)=zinds(l,c);
      int q2=vq(i)=zindq(l,c);
      double angrad2=anf.atf(l,c).rotmap-tef(s2,q2).rotmap; // em radianos
      match(i)=mncc(an,te(s2,q2),l,c,angrad2);
      vang(i)=angrad2;
    }
    mordena3(match,v,vang,vesc,vq,0,match.n()-1); 
    p=v[v.size()-1];
    angrad=-vang(v.size()-1);
    s=vesc(v.size()-1);
    q=vq(v.size()-1);
  } else {
    p=argmin(zbuf);
    int l=p.l(); int c=p.c();
    s=zinds(l,c);
    q=zindq(l,c);
    angrad=-anf.atf(l,c).rotmap+tef(s).rotmap; // em radianos
  }
  int tnl2=arredonda(escinic*pow(passoesc,s)*tnl);
  relogioretangular2(sai,p.l(),p.c(),angrad,tnl2,tnl2,q,s); // No futuro, trocar s por q
  pvImp(sai,argv[2]);

  int t3=centseg();
  printf("Conv_fea_a=%d procura=%d total=%d\n",t2-t1,t3-t2,t3-t1);
}

void msvhfi1(int argc, char** argv)
{
  if (argc<15) {
    printf("MSVHfi1: Scale invariant VHough. Cd Q aparece 1 vez em A\n");
    printf("SVHfi1 A.tga P.tga nrad ncir prec detec nte ncand escinic escfim nesc tnl separacao Q*.tga\n");
    printf("  Ex: a.tga p.tga 4 4 0.1 1.6 4 10 0.7 1.4 8 31 11 q*.tga\n");
    erro("Error: Invalid number of arguments");
  }

  int t1=centseg();
  IMGGRY an; pvLe(an,argv[1]);

  int nrad;
  if (sscanf(argv[3],"%d",&nrad)!=1) erro("Erro: Leitura nrad");
  if (nrad<1) erro("Erro: nrad<1");

  int ncir;
  if (sscanf(argv[4],"%d",&ncir)!=1) erro("Erro: Leitura ncir");
  if (ncir<0) erro("Erro: ncir<0");

  double precisao;
  if (sscanf(argv[5],"%lf",&precisao)!=1) erro("Erro: Leitura precisao");
  if (precisao<0.0 || 1.0<precisao) erro("Erro: Precisao fora de [0,1]");

  double detec;
  if (sscanf(argv[6],"%lf",&detec)!=1) erro("Erro: Leitura detec");
  if (detec<1.0) erro("Erro: detec<1.0");

  int nte;
  if (sscanf(argv[7],"%d",&nte)!=1) erro("Erro: Leitura nte");
  if (nte<0) erro("Erro: nte<0");

  int ncand;
  if (sscanf(argv[8],"%d",&ncand)!=1) erro("Erro: Leitura ncand");
  if (ncand<0) erro("Erro: ncand<0");

  double escinic;
  if (sscanf(argv[9],"%lf",&escinic)!=1) erro("Erro: Leitura escinic");
  if (escinic<0.0) erro("Erro: escinic<0.0");

  double escfim;
  if (sscanf(argv[10],"%lf",&escfim)!=1) erro("Erro: Leitura escfim");
  if (escfim<escinic) erro("Erro: escfim<escinic");

  int nesc;
  if (sscanf(argv[11],"%d",&nesc)!=1) erro("Erro: Leitura nesc");
  if (nesc<1) erro("Erro: nesc<1");
  double passoesc=exp(log(escfim/escinic)/(nesc-1));

  int tnl;
  if (sscanf(argv[12],"%d",&tnl)!=1) erro("Erro: Leitura tnl");
  if (tnl%2==0 || tnl<=0) erro("Erro: tnl deve ser impar");
  int raio=tnl/2;

  int separacao;
  if (sscanf(argv[13],"%d",&separacao)!=1) erro("Erro: Leitura separacao");
  if (separacao<0 || tnl<separacao) erro("Erro: separacao<0 || tnl<separacao");

  IMGGRY anv(an.nl()-tnl+1,an.nc()-tnl+1);
  for (int l=0; l<anv.nl(); l++)
    for (int c=0; c<anv.nc(); c++) {
      anv(l,c)=an(l+raio,c+raio);
  }

  double deteccao=detec*precisao;
  PESOTOTAL pesototal(nrad,ncir);

  I3DCPX qux; // nrad+ncir, tnl, tnl
  geraqux(qux,nrad,ncir,tnl);

  I3DCPX anc;
  convqux(an,qux,anc);

  MFEA anf; // 2*nrad+2*ncir, quv.nl(), quv.nc()
  IMGGRY valido;
  gerafea(anc,anf,valido,nrad,ncir,tnl);

  IMGGRY qu;
  IMGCOR sai=an;
  int arginic=14;

  int t2=centseg();
  for (int q=arginic; q<argc; q++) {
    pvLe(qu,argv[q]);
    printf("Processando %s\n",argv[q]);

    // Geracao dos Qs e templates
    VETOR<IMGGRY> qesc(nesc);
    for (int s=0; s<nesc; s++) {
      double esc=escinic*pow(passoesc,s);
      int nl=arredonda(esc*qu.nl());
      int nc=arredonda(esc*qu.nc());
      qesc(s)=zoomnp(qu,nl,nc);
    }

    VETOR<PONTOI2> posesc(nesc);
    VETOR<double>  houesc(nesc);
    VETOR<double>  anaesc(nesc);

    for (int s=0; s<nesc; s++) {
      IMGGRY& qus=qesc(s);

      IMGGRY quv(qus.nl()-tnl+1,qus.nc()-tnl+1);
      for (int l=0; l<quv.nl(); l++)
        for (int c=0; c<quv.nc(); c++)
          quv(l,c)=qus(l+raio,c+raio);

      I3DCPX quc;
      convqux(qus,qux,quc);

      MFEA quf; // 2*nrad+2*ncir, quv.nl(), quv.nc()
      gerafea(quc,quf,valido,nrad,ncir,tnl);

      stablecand(quf,valido,precisao,pesototal);

      IMGFLT f=distgeo(valido);
      //pvMostra(IMGGRY(normaliza(f)),"valido distgeo",1);
      vector<PONTOI2> vt=kmax(f,nte,separacao);
      if (vt.size()==0) { // erro("Erro: Sem ponto estavel");
        printf("Warning: Escala=%d. #templates procurados=%d nenhum template achado.\n",s,nte);
        continue;
      } else if (nte!=int(vt.size())) printf("Warning: Escala=%d #templates procurados=%d achados=%d\n",s,nte,vt.size());

      IMGFLT mat;
      FEA tef;

      MATRIZ<PONTOI2> pos(nte, ncand,PONTOI2(maxint,maxint));
      IMGCPX  ana(nte, ncand,CPX(0,0));
      IMGFLT  hou(nte, ncand,-infinito);

      double passo=max(tnl/5,3); // tolerancia
      for (unsigned i=0; i<vt.size(); i++) {
        tef=quf(vt[i].l(),vt[i].c());
        int dl=quv.lc()-vt[i].l(); int dc=quv.cc()-vt[i].c();
        //procura1(anf,tef,mat,deteccao,pesototal);
        // Acho que nao pode usar procura1 aqui, pois vai procurar ncand matchings
        procura(anf,tef,mat,pesototal); 

        vector<PONTOI2> vh=kargmin(mat,ncand);
        if (ncand!=int(vh.size())) printf("Warning: Escala=%d template=%d #cand procurados=%d achados=%d\n",s,i,ncand,vh.size());

        for (unsigned j=0; j<vh.size(); j++) {
          int l=vh[j].l(); int c=vh[j].c();
          double peso=1.0-mat.atf(l,c);
          double angulo=-anf.atf(l,c).rotmap+tef.rotmap;
          double co=cos(angulo); double si=sin(angulo);
          int ndl=arredonda( co*dl+si*dc);
          int ndc=arredonda(-si*dl+co*dc);
          hou.atf(i,j)=peso;
          ana.atf(i,j)=CPX(peso*co,peso*si);
          pos.atf(i,j)=PONTOI2(l+ndl,c+ndc);
        }
      }

      double maxsomapeso=-infinito;
      PONTO2 maxmedia(0.0,0.0);
      CPX maxana(0.0,0.0);
      int maxconta=0;
      for (int c=0; c<pos.nc(); c++) {
        for (int l=0; l<pos.nl(); l++) { // Tem que fazer loop c/l para pegar melhor matching
          double somapeso=hou.atf(l,c);
          PONTO2 somaponto=somapeso*PONTO2(pos.atf(l,c));
          CPX somaana=ana.atf(l,c);
          int conta=1;
          for (int l2=0; l2<pos.nl(); l2++) {
            if (l==l2) break;
            for (int c2=0; c2<pos.nc(); c2++) {
              if (distancia(pos.atf(l,c),pos.atf(l2,c2))<passo) {
                somapeso += hou.atf(l2,c2);
                somaponto = somaponto + hou.atf(l2,c2)*PONTO2(pos.atf(l2,c2));
                somaana += ana.atf(l2,c2); // Nao esta multiplicando peso
                conta++;
                break;
              }
            }
          }
          if (somapeso>maxsomapeso) {
            maxsomapeso=somapeso;
            maxmedia=(1.0/somapeso)*somaponto;
            maxana=somaana;
            maxconta=conta;
            if (maxconta==pos.nl()) goto saida;
          }
        }
      }
      saida:

      anaesc(s)=arg(maxana);
      int l=arredonda(maxmedia.l())+raio; int c=arredonda(maxmedia.c())+raio; posesc(s)=PONTOI2(l,c);
      houesc(s)=maxsomapeso;
    }
    int maxs=argmax(houesc);
    relogioretangular2(sai,posesc(maxs).l(),posesc(maxs).c(),anaesc(maxs),qesc(maxs).nl(),qesc(maxs).nc(),q-arginic+1);
    //relogioretangular2(sai,posesc(maxs).l()+raio,posesc(maxs).c()+raio,anaesc(maxs),qesc(maxs).nl(),qesc(maxs).nc(),maxs);
  }
  pvImp(sai,argv[2]);

  int t3=centseg();
  printf("pre=%d busca=%d total=%d\n",t2-t1,t3-t2,t3-t1);
}

void afvhfi1(int argc, char** argv)
{
  if (argc<15) {
    printf("AfVHfi1: Affine invariant VHough. Cd Q aparece 1 vez em A\n");
    printf("  Deve fornecer Qs distorcidos com fundo branco (255)\n");
    printf("AfVHfi1 A.tga P.tga nrad ncir prec detec nte ncand escinic escfim nesc tnl separacao Q*.tga\n");
    printf("  Ex: a.tga p.tga 4 4 0.1 1.6 4 10 0.7 1.4 4 31 11 q*.tga\n");
    erro("Error: Invalid number of arguments");
  }

  int t1=centseg();
  IMGGRY an; pvLe(an,argv[1]);

  int nrad;
  if (sscanf(argv[3],"%d",&nrad)!=1) erro("Erro: Leitura nrad");
  if (nrad<1) erro("Erro: nrad<1");

  int ncir;
  if (sscanf(argv[4],"%d",&ncir)!=1) erro("Erro: Leitura ncir");
  if (ncir<0) erro("Erro: ncir<0");

  double precisao;
  if (sscanf(argv[5],"%lf",&precisao)!=1) erro("Erro: Leitura precisao");
  if (precisao<0.0 || 1.0<precisao) erro("Erro: Precisao fora de [0,1]");

  double detec;
  if (sscanf(argv[6],"%lf",&detec)!=1) erro("Erro: Leitura detec");
  if (detec<1.0) erro("Erro: detec<1.0");

  int nte;
  if (sscanf(argv[7],"%d",&nte)!=1) erro("Erro: Leitura nte");
  if (nte<0) erro("Erro: nte<0");

  int ncand;
  if (sscanf(argv[8],"%d",&ncand)!=1) erro("Erro: Leitura ncand");
  if (ncand<0) erro("Erro: ncand<0");

  double escinic;
  if (sscanf(argv[9],"%lf",&escinic)!=1) erro("Erro: Leitura escinic");
  if (escinic<0.0) erro("Erro: escinic<0.0");

  double escfim;
  if (sscanf(argv[10],"%lf",&escfim)!=1) erro("Erro: Leitura escfim");
  if (escfim<escinic) erro("Erro: escfim<escinic");

  int nesc;
  if (sscanf(argv[11],"%d",&nesc)!=1) erro("Erro: Leitura nesc");
  if (nesc<1) erro("Erro: nesc<1");
  double passoesc=exp(log(escfim/escinic)/(nesc-1));

  int tnl;
  if (sscanf(argv[12],"%d",&tnl)!=1) erro("Erro: Leitura tnl");
  if (tnl%2==0 || tnl<=0) erro("Erro: tnl deve ser impar");
  int raio=tnl/2;

  int separacao;
  if (sscanf(argv[13],"%d",&separacao)!=1) erro("Erro: Leitura separacao");
  if (separacao<0 || tnl<separacao) erro("Erro: separacao<0 || tnl<separacao");

  IMGGRY anv(an.nl()-tnl+1,an.nc()-tnl+1);
  for (int l=0; l<anv.nl(); l++)
    for (int c=0; c<anv.nc(); c++) {
      anv(l,c)=an(l+raio,c+raio);
  }

  double deteccao=detec*precisao;
  PESOTOTAL pesototal(nrad,ncir);

  I3DCPX qux; // nrad+ncir, tnl, tnl
  geraqux(qux,nrad,ncir,tnl);

  I3DCPX anc;
  convqux(an,qux,anc);

  MFEA anf; // 2*nrad+2*ncir, quv.nl(), quv.nc()
  IMGGRY valido;
  gerafea(anc,anf,valido,nrad,ncir,tnl);

  IMGGRY qu;
  IMGCOR sai=an;
  int arginic=14;

  int maxindq=0;
  int maxscaq=0;
  double maxhouq=-infinito;
  double maxanaq=0.0;

  int t2=centseg();
  for (int q=arginic; q<argc; q++) {
    pvLe(qu,argv[q]);
    printf("Processando %s\n",argv[q]);

    // Geracao dos Qs e templates
    VETOR<IMGGRY> qesc(nesc);
    for (int s=0; s<nesc; s++) {
      double esc=escinic*pow(passoesc,s);
      int nl=arredonda(esc*qu.nl());
      int nc=arredonda(esc*qu.nc());
      qesc(s)=zoomnp(qu,nl,nc);
    }

    VETOR<PONTOI2> posesc(nesc); 
    VETOR<double>  houesc(nesc); 
    VETOR<double>  anaesc(nesc); 

    for (int s=0; s<nesc; s++) {
      IMGGRY& qus=qesc(s);

      IMGGRY quv(qus.nl()-tnl+1,qus.nc()-tnl+1);
      for (int l=0; l<quv.nl(); l++)
        for (int c=0; c<quv.nc(); c++)
          quv(l,c)=qus(l+raio,c+raio);

      I3DCPX quc;
      convqux(qus,qux,quc);

      MFEA quf; // 2*nrad+2*ncir, quv.nl(), quv.nc()
      gerafea(quc,quf,valido,nrad,ncir,tnl);
      if (ignora255) Ignora255(qus,valido,tnl);

      stablecand(quf,valido,precisao,pesototal);

      IMGFLT f=distgeo(valido);
      //pvMostra(IMGGRY(normaliza(f)),"valido distgeo",1);
      vector<PONTOI2> vt=kmax(f,nte,separacao);
      if (vt.size()==0) { // erro("Erro: Sem ponto estavel");
        printf("Warning: Escala=%d. #templates procurados=%d nenhum template achado.\n",s,nte);
        continue;
      } else if (nte!=int(vt.size())) printf("Warning: Escala=%d #templates procurados=%d achados=%d\n",s,nte,vt.size());

      IMGFLT mat;
      FEA tef;

      MATRIZ<PONTOI2> pos(nte, ncand,PONTOI2(maxint,maxint)); 
      IMGCPX  ana(nte, ncand,CPX(0,0));
      IMGFLT  hou(nte, ncand,-infinito);

      double passo=max(tnl/5,3); // tolerancia
      for (unsigned i=0; i<vt.size(); i++) {
        tef=quf(vt[i].l(),vt[i].c());
        int dl=quv.lc()-vt[i].l(); int dc=quv.cc()-vt[i].c();
        //procura1(anf,tef,mat,deteccao,pesototal); 
        // Acho que nao pode usar procura1 aqui, pois vai procurar ncand matchings
        procura(anf,tef,mat,pesototal); 

        vector<PONTOI2> vh=kargmin(mat,ncand);
        if (ncand!=int(vh.size())) 
           printf("Warning: Escala=%d template=%d #cand procurados=%d achados=%d\n",s,i,ncand,vh.size());

        for (unsigned j=0; j<vh.size(); j++) {
          int l=vh[j].l(); int c=vh[j].c();
          double peso=1.0-mat.atf(l,c);
          double angulo=-anf.atf(l,c).rotmap+tef.rotmap;
          double co=cos(angulo); double si=sin(angulo);
          int ndl=arredonda( co*dl+si*dc);
          int ndc=arredonda(-si*dl+co*dc);
          hou.atf(i,j)=peso;
          ana.atf(i,j)=CPX(peso*co,peso*si);
          pos.atf(i,j)=PONTOI2(l+ndl,c+ndc);
        }
      }

      double maxsomapeso=-infinito;
      PONTO2 maxmedia(0.0,0.0);
      CPX maxana(0.0,0.0);
      int maxconta=0;
      for (int c=0; c<pos.nc(); c++) {
        for (int l=0; l<pos.nl(); l++) { // Tem que fazer loop c/l para pegar melhor matching
          double somapeso=hou.atf(l,c);
          PONTO2 somaponto=somapeso*PONTO2(pos.atf(l,c));
          CPX somaana=ana.atf(l,c);
          int conta=1;
          for (int l2=0; l2<pos.nl(); l2++) {
            if (l==l2) break;
            for (int c2=0; c2<pos.nc(); c2++) {
              if (distancia(pos.atf(l,c),pos.atf(l2,c2))<passo) {
                somapeso += hou.atf(l2,c2);
                somaponto = somaponto + hou.atf(l2,c2)*PONTO2(pos.atf(l2,c2));
                somaana += ana.atf(l2,c2); // Nao esta multiplicando peso
                conta++;
                break;
              }
            }
          }
          if (somapeso>maxsomapeso) {
            maxsomapeso=somapeso;
            maxmedia=(1.0/somapeso)*somaponto;
            maxana=somaana;
            maxconta=conta;
            if (maxconta==pos.nl()) goto saida;
          }
        }
      }
      saida:

      anaesc(s)=arg(maxana);
      int l=arredonda(maxmedia.l())+raio; int c=arredonda(maxmedia.c())+raio; posesc(s)=PONTOI2(l,c);
      houesc(s)=maxsomapeso;
    }
    int maxs=argmax(houesc); // houesc indica a qualidade do matching.
    double maxh=houesc(maxs);
    if (maxhouq<maxh) { maxhouq=maxh; maxindq=q-arginic; maxscaq=maxs; maxanaq=anaesc(maxs); }

    relogioretangular2(sai,posesc(maxs).l(),posesc(maxs).c(),anaesc(maxs),qesc(maxs).nl(),qesc(maxs).nc(),q-arginic);
    printf("#arq=%d qualidade=%f maxscale=%d(%f) nomearq=%s\n",q-arginic,maxh,maxs,escinic*pow(passoesc,maxs),argv[q]);
  }
  pvImp(sai,argv[2]);
  printf("Melhor matching indice=%d qualidade=%f maxscale=%d(%f) ang=%d nomearq=%s\n",
         maxindq,maxhouq,maxscaq,escinic*pow(passoesc,maxscaq),arredonda(rad2deg(maxanaq)),argv[maxindq+arginic]);
  int t3=centseg();
  printf("pre=%d busca=%d total=%d\n",t2-t1,t3-t2,t3-t1);
}

IMGGRY CombineVertically(IMGGRY im1, IMGGRY im2)
{ int rows, cols, r, c;
  IMGGRY result;

  rows = im1.nl() + im2.nl();
  cols = max(im1.nc(), im2.nc());
  result.resize(rows, cols);
  result.fill(255);

  for (r = 0; r < im1.nl(); r++)
    for (c = 0; c < im1.nc(); c++)
      result(r,c) = im1(r,c);
  for (r = 0; r < im2.nl(); r++)
    for (c = 0; c < im2.nc(); c++)
      result(r + im1.nl(),c) = im2(r,c);

  return result;
}

void vhfi1(int argc, char** argv)
{ if (argc!=12) {
    printf("VHfi1: Forapro VHough com linhas entre A e Q. Cd Q aparece 1 vez em A\n");
    printf("  Deve fornecer Q com fundo branco (255)\n");
    printf("VHfi1 A.tga P.tga nrad ncir prec nte ncand esc tnl separacao Q.tga\n");
    printf("  Ex: a.tga p.tga 4    4    0.1  4   10    1.0 31  11        q.tga\n");
    erro("Error: Invalid number of arguments");
  }

  IMGGRY an; pvLe(an,argv[1]);

  int nrad;
  if (sscanf(argv[3],"%d",&nrad)!=1) erro("Erro: Leitura nrad");
  if (nrad<1) erro("Erro: nrad<1");

  int ncir;
  if (sscanf(argv[4],"%d",&ncir)!=1) erro("Erro: Leitura ncir");
  if (ncir<0) erro("Erro: ncir<0");

  double precisao;
  if (sscanf(argv[5],"%lf",&precisao)!=1) erro("Erro: Leitura precisao");
  if (precisao<0.0 || 1.0<precisao) erro("Erro: Precisao fora de [0,1]");

  int nte;
  if (sscanf(argv[6],"%d",&nte)!=1) erro("Erro: Leitura nte");
  if (nte<0) erro("Erro: nte<0");

  int ncand;
  if (sscanf(argv[7],"%d",&ncand)!=1) erro("Erro: Leitura ncand");
  if (ncand<0) erro("Erro: ncand<0");

  double esc;
  if (sscanf(argv[8],"%lf",&esc)!=1) erro("Erro: Leitura esc");
  if (esc<0.0) erro("Erro: esc<0.0");

  int tnl;
  if (sscanf(argv[9],"%d",&tnl)!=1) erro("Erro: Leitura tnl");
  if (tnl%2==0 || tnl<=0) erro("Erro: tnl deve ser impar");
  int raio=tnl/2;

  int separacao;
  if (sscanf(argv[10],"%d",&separacao)!=1) erro("Erro: Leitura separacao");
  if (separacao<0 || tnl<separacao) erro("Erro: separacao<0 || tnl<separacao");

  IMGGRY anv(an.nl()-tnl+1,an.nc()-tnl+1);
  for (int l=0; l<anv.nl(); l++)
    for (int c=0; c<anv.nc(); c++) {
      anv(l,c)=an(l+raio,c+raio);
  }

  PESOTOTAL pesototal(nrad,ncir);

  I3DCPX qux; // nrad+ncir, tnl, tnl
  geraqux(qux,nrad,ncir,tnl);

  I3DCPX anc;
  convqux(an,qux,anc);

  MFEA anf; // 2*nrad+2*ncir, quv.nl(), quv.nc()
  IMGGRY valido;
  gerafea(anc,anf,valido,nrad,ncir,tnl);

  IMGGRY qu;

  pvLe(qu,argv[11]);

  // Geracao dos Qs e templates
  int nl=arredonda(esc*qu.nl());
  int nc=arredonda(esc*qu.nc());
  IMGGRY qus=zoomnp(qu,nl,nc);

  PONTOI2 posesc; 
  double  houesc; 
  double  anaesc; 

  IMGGRY quv(qus.nl()-tnl+1,qus.nc()-tnl+1);
  for (int l=0; l<quv.nl(); l++)
    for (int c=0; c<quv.nc(); c++)
      quv(l,c)=qus(l+raio,c+raio);

  I3DCPX quc;
  convqux(qus,qux,quc);

  MFEA quf; // 2*nrad+2*ncir, quv.nl(), quv.nc()
  gerafea(quc,quf,valido,nrad,ncir,tnl);
  if (ignora255) Ignora255(qus,valido,tnl);

  stablecand(quf,valido,precisao,pesototal);

  IMGFLT f=distgeo(valido);
  //pvMostra(IMGGRY(normaliza(f)),"valido distgeo",1);
  vector<PONTOI2> vt=kmax(f,nte,separacao);
  if (vt.size()==0) { // erro("Erro: Sem ponto estavel");
    printf("Erro: #templates procurados=%d nenhum template achado.\n",nte);
    exit(0);
  } else if (nte!=int(vt.size())) 
    printf("Warning: #templates procurados=%d achados=%d\n",nte,vt.size());

  IMGFLT mat;
  FEA tef;

  MATRIZ<PONTOI2> pos(nte, ncand,PONTOI2(maxint,maxint)); 
  MATRIZ<PONTOI2> posori(nte, ncand,PONTOI2(maxint,maxint)); 
  IMGCPX  ana(nte, ncand,CPX(0,0));
  IMGFLT  hou(nte, ncand,-infinito);

  double passo=max(tnl/5,3); // tolerancia
  for (unsigned i=0; i<vt.size(); i++) {
    tef=quf(vt[i].l(),vt[i].c());
    int dl=quv.lc()-vt[i].l(); int dc=quv.cc()-vt[i].c();
    //procura1(anf,tef,mat,deteccao,pesototal); 
    // Acho que nao pode usar procura1 aqui, pois vai procurar ncand matchings
    procura(anf,tef,mat,pesototal); 

    vector<PONTOI2> vh=kargmin(mat,ncand);
    if (ncand!=int(vh.size())) 
       printf("Warning: Template=%d #cand procurados=%d achados=%d\n",i,ncand,vh.size());

    for (unsigned j=0; j<vh.size(); j++) {
      int l=vh[j].l(); int c=vh[j].c();
      posori.atf(i,j)=PONTOI2(l,c);
      double peso=1.0-mat.atf(l,c);
      double angulo=-anf.atf(l,c).rotmap+tef.rotmap;
      double co=cos(angulo); double si=sin(angulo);
      int ndl=arredonda( co*dl+si*dc);
      int ndc=arredonda(-si*dl+co*dc);
      hou.atf(i,j)=peso;
      ana.atf(i,j)=CPX(peso*co,peso*si);
      pos.atf(i,j)=PONTOI2(l+ndl,c+ndc);
    }
  }

  double maxsomapeso=-infinito;
  PONTO2 maxmedia(0.0,0.0);
  CPX maxana(0.0,0.0);
  int maxconta=0;
  VETOR<int> maxind(pos.nl(),-1);
  VETOR<int> ind(pos.nl());
  for (int c=0; c<pos.nc(); c++) {
    ind.fill(-1);
    for (int l=0; l<pos.nl(); l++) { // Tem que fazer loop c/l para pegar melhor matching
      ind(l)=c;
      double somapeso=hou.atf(l,c);
      PONTO2 somaponto=somapeso*PONTO2(pos.atf(l,c));
      CPX somaana=ana.atf(l,c);
      int conta=1;
      for (int l2=0; l2<pos.nl(); l2++) {
        if (l==l2) break;
        for (int c2=0; c2<pos.nc(); c2++) {
          if (distancia(pos.atf(l,c),pos.atf(l2,c2))<passo) {
            somapeso += hou.atf(l2,c2);
            somaponto = somaponto + hou.atf(l2,c2)*PONTO2(pos.atf(l2,c2));
            somaana += ana.atf(l2,c2); // Nao esta multiplicando peso
            conta++;
            ind(l2)=c2;
            break;
          }
        }
      }
      if (somapeso>maxsomapeso) {
        maxsomapeso=somapeso;
        maxmedia=(1.0/somapeso)*somaponto;
        maxana=somaana;
        maxconta=conta;
        maxind=ind;
        if (maxconta==pos.nl()) goto saida;
      }
    }
  }
  saida:

  anaesc=arg(maxana);
  int l=arredonda(maxmedia.l())+raio; int c=arredonda(maxmedia.c())+raio; 

  IMGCOR sai=CombineVertically(qus,an);
  for (unsigned i=0; i<vt.size(); i++) {
    int li=vt[i].l()+raio; int ci=vt[i].c()+raio;
    int j=maxind(i);
    if (j<0) continue;
    int lf=posori(i,j).l()+qus.nl()+raio; int cf=posori(i,j).c()+raio;
    reta(sai,li,ci,lf,cf,COR(255,0,0));
    ponto(sai,li,ci,3,COR(0,255,0));
    ponto(sai,lf,cf,3,COR(0,255,0));
  }
  relogioretangular2(sai,qus.lc(),qus.cc(),0,qus.nl()/3,qus.nc()/3);
  relogioretangular2(sai,l+qus.nl(),c,anaesc,qus.nl()/3,qus.nc()/3);
  //pvMostra(sai);
  pvImp(sai,argv[2]);

  printf("qualidade=%f l=%d c=%d ang=%d escala=%f\n",maxsomapeso,l,c,arredonda(rad2deg(anaesc)),esc);

}

//<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< main <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
const char about[]=
"< Forapro: Forapro Template Matching that ignores white pixels, version 3.73>\n"
"_______________________________________________________________________________\n";

const char help[]=
"Programas:\n"
"GeraTe  -Gera n templates estaveis a partir de query\n"
"TeRuim  -Gera um template propositalmente ruim a partir de query\n"
"Features-Gera rotation-discriminating and -invariant features.\n"
"FeaDist -Matching AxT com features pre-calculados.\n"
"MatDist -(TeMatch) Matching AxT. Saida grayscale matching distance.\n"
"MTefi1  -Matching A x varios Ts com NCC. A contem 1 T. Nao calcula estavel\n"
"MVHfi1  -Matching A x varios Qs com VHough. Gera Ts. Cd Q aparece 1 vez em A\n"
"MSTefi1 -Scale invariant MTefi com NCC. Gera Ts. Cd Q aparece 1 vez em A\n"
"AfTefi1 -Affine invariant MTefi com NCC. Gera Ts. Cd Q aparece 1 vez em A\n"
"MSVHfi1 -Scale invariant VHough. Gera Ts. Cd Q aparece 1 vez em A\n"
"AfVHfi1 -Affine invariant VHough. Gera Ts. Cd Q aparece 1 vez em A\n"
"VHfi1   -VHough com linhas entre A e Q. Gera Ts. Cd Q aparece 1 vez em A\n"
"...............................................................................\n";

int main(int argc, char** argv)
{
  printf(about);
  if (argc<2) {
    cout << help;
    erro("Erro: Numero de argumentos invalido");
  } else {
    string comando=strlwr(argv[1]);
    if      (comando=="gerate")  gerate(argc-1,&argv[1]);
    else if (comando=="teruim")  teruim(argc-1,&argv[1]);
    else if (comando=="feadist") feadist(argc-1,&argv[1]);
    else if (comando=="features") features(argc-1,&argv[1]);
    else if (comando=="matdist") matdist(argc-1,&argv[1]);
    else if (comando=="mtefi1")  mtefi1(argc-1,&argv[1]);
    else if (comando=="mvhfi1")  mvhfi1(argc-1,&argv[1]);
    else if (comando=="mstefi1") mstefi1(argc-1,&argv[1]);
    else if (comando=="aftefi1") aftefi1(argc-1,&argv[1]);
    else if (comando=="msvhfi1") msvhfi1(argc-1,&argv[1]);
    else if (comando=="afvhfi1") afvhfi1(argc-1,&argv[1]);
    else if (comando=="vhfi1")   vhfi1(argc-1,&argv[1]);
    else erro("Erro: Programa inexistente ",comando.c_str());
  }
}