/***************************************************************************
Haar 2d transform implemented in C/C++ to speed things up
-------------------
Implemented for imgSeek by nieder|at|mail.ru
based off of
***************************************************************************
* Wavelet algorithms, metric and query ideas based on the paper *
* Fast Multiresolution Image Querying *
* by Charles E. Jacobs, Adam Finkelstein and David H. Salesin. *
* <http://www.cs.washington.edu/homes/salesin/abstracts.html> *
***************************************************************************
Version from imgSeek Copyright (C) 2003 Ricardo Niederberger Cabral
XS version derived & Copyright (C) 2005 Simon Cozens
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 "EXTERN.h"
#include "perl.h"
#include <queue>
#include "haar.h"
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
inline double *absarray(double* a) {
double * b;
New(200, b, 16384, double);
int i;
for (i=0;i<16384;i++) b[i] = fabs(a[i]);
return b;
}
inline void truncq(double* a,double* b,double limit,int* data) {
/* set every element on data to 0 (if < limit) or to a positive/negative int when abs(value) > limit.
a is the absolute value of array b
*/
int i,c=0;
Zero(data,40, int);
for (i=0;i<16384;i++) {
if (a[i]>limit) {
if (b[i]>0) data[c]=i;
else data[c]=-1*i;
c++;
if (c==40) { return; }
}
}
}
void transformChar(unsigned char* c1,unsigned char* c2,unsigned char* c3,double* a,double* b,double* c) {
// do the Haar tensorial 2d transform itself
// here input RGB data is on unsigned char arrays
int i,h,j,k;
double *d1, *d2, *d3, *Ab1, *Ab2, *Ab3;
New(200, d1, 16384, double);
New(200, d2, 16384, double);
New(200, d3, 16384, double);
New(200, Ab1, 128, double);
New(200, Ab2, 128, double);
New(200, Ab3, 128, double);
/* RGB -> YIQ */
for (i=0;i<16384;i++) {
d1[i]=(c1[i]*0.299+c2[i]*0.587+c3[i]*0.114)/256.0;
d2[i]=(c1[i]*0.596+c2[i]*(-0.274)+c3[i]*(-0.322))/256.0;
d3[i]=(c1[i]*0.212+c2[i]*(-0.523)+c3[i]*0.311)/256.0;
}
// decompose rows
for (i=0;i<128;i++) {
h=128;
for (j=0;j<128;j++) {
d1[i*128+j] /= 11.314;
d2[i*128+j] /= 11.314;
d3[i*128+j] /= 11.314;
}
while(h>1) {
h/=2;
for (k=0;k<h;k++) {
Ab1[k]=(d1[i*128+2*k]+d1[i*128+2*k+1])/1.414;
Ab2[k]=(d2[i*128+2*k]+d2[i*128+2*k+1])/1.414;
Ab3[k]=(d3[i*128+2*k]+d3[i*128+2*k+1])/1.414;
Ab1[k+h]=(d1[i*128+2*k]-d1[i*128+2*k+1])/1.414;
Ab2[k+h]=(d2[i*128+2*k]-d2[i*128+2*k+1])/1.414;
Ab3[k+h]=(d3[i*128+2*k]-d3[i*128+2*k+1])/1.414;
}
memcpy(d1+i*128,Ab1,2*h*sizeof(double));
memcpy(d2+i*128,Ab2,2*h*sizeof(double));
memcpy(d3+i*128,Ab3,2*h*sizeof(double));
}
}
// decompose cols
for (i=0;i<128;i++) {
h=128;
for (j=0;j<128;j++) {
d1[j*128+i] /= 11.314;
d2[j*128+i] /= 11.314;
d3[j*128+i] /= 11.314;
}
while(h>1) {
h/=2;
for (k=0;k<h;k++) {
Ab1[k]=(d1[i+2*k*128]+d1[i+(2*k+1)*128])/1.414;
Ab2[k]=(d2[i+2*k*128]+d2[i+(2*k+1)*128])/1.414;
Ab3[k]=(d3[i+2*k*128]+d3[i+(2*k+1)*128])/1.414;
Ab1[k+h]=(d1[i+2*k*128]-d1[i+(2*k+1)*128])/1.414;
Ab2[k+h]=(d2[i+2*k*128]-d2[i+(2*k+1)*128])/1.414;
Ab3[k+h]=(d3[i+2*k*128]-d3[i+(2*k+1)*128])/1.414;
}
for (j=0;j<2*h;j++) {
d1[j*128+i]=Ab1[j];
d2[j*128+i]=Ab2[j];
d3[j*128+i]=Ab3[j];
}
}
}
memcpy(a,d1,16384*sizeof(double));
memcpy(b,d2,16384*sizeof(double));
memcpy(c,d3,16384*sizeof(double));
Safefree(d1);Safefree(d2);Safefree(d3);
Safefree(Ab1);Safefree(Ab2);Safefree(Ab3);
}
void transform(double* a,double* b,double* c) {
// do the Haar tensorial 2d transform itself
// here input RGB data is on double char arrays
int i,h,j,k;
double *d1, *d2, *d3, *Ab1, *Ab2, *Ab3;
New(200, d1, 16384, double);
New(200, d2, 16384, double);
New(200, d3, 16384, double);
New(200, Ab1, 128, double);
New(200, Ab2, 128, double);
New(200, Ab3, 128, double);
/* RGB -> YIQ */
for (i=0;i<16384;i++) {
d1[i]=(a[i]*0.299+b[i]*0.587+c[i]*0.114)/256.0;
d2[i]=(a[i]*0.596+b[i]*(-0.274)+c[i]*(-0.322))/256.0;
d3[i]=(a[i]*0.212+b[i]*(-0.523)+c[i]*0.311)/256.0;
}
// decompose rows
for (i=0;i<128;i++) {
h=128;
for (j=0;j<128;j++) {
d1[i*128+j] /= 11.314;
d2[i*128+j] /= 11.314;
d3[i*128+j] /= 11.314;
}
while(h>1) {
h/=2;
for (k=0;k<h;k++) {
Ab1[k]=(d1[i*128+2*k]+d1[i*128+2*k+1])/1.414;
Ab2[k]=(d2[i*128+2*k]+d2[i*128+2*k+1])/1.414;
Ab3[k]=(d3[i*128+2*k]+d3[i*128+2*k+1])/1.414;
Ab1[k+h]=(d1[i*128+2*k]-d1[i*128+2*k+1])/1.414;
Ab2[k+h]=(d2[i*128+2*k]-d2[i*128+2*k+1])/1.414;
Ab3[k+h]=(d3[i*128+2*k]-d3[i*128+2*k+1])/1.414;
}
memcpy(d1+i*128,Ab1,2*h*sizeof(double));
memcpy(d2+i*128,Ab2,2*h*sizeof(double));
memcpy(d3+i*128,Ab3,2*h*sizeof(double));
}
}
// decompose cols
for (i=0;i<128;i++) {
h=128;
for (j=0;j<128;j++) {
d1[j*128+i] /= 11.314;
d2[j*128+i] /= 11.314;
d3[j*128+i] /= 11.314;
}
while(h>1) {
h/=2;
for (k=0;k<h;k++) {
Ab1[k]=(d1[i+2*k*128]+d1[i+(2*k+1)*128])/1.414;
Ab2[k]=(d2[i+2*k*128]+d2[i+(2*k+1)*128])/1.414;
Ab3[k]=(d3[i+2*k*128]+d3[i+(2*k+1)*128])/1.414;
Ab1[k+h]=(d1[i+2*k*128]-d1[i+(2*k+1)*128])/1.414;
Ab2[k+h]=(d2[i+2*k*128]-d2[i+(2*k+1)*128])/1.414;
Ab3[k+h]=(d3[i+2*k*128]-d3[i+(2*k+1)*128])/1.414;
}
for (j=0;j<2*h;j++) {
d1[j*128+i]=Ab1[j];
d2[j*128+i]=Ab2[j];
d3[j*128+i]=Ab3[j];
}
}
}
memcpy(a,d1,16384*sizeof(double));
memcpy(b,d2,16384*sizeof(double));
memcpy(c,d3,16384*sizeof(double));
Safefree(d1);Safefree(d2);Safefree(d3);Safefree(Ab1);Safefree(Ab2);Safefree(Ab3);
}
int calcHaar(double* cdata1,double* cdata2,double* cdata3,int* sig1,int* sig2,int* sig3,double * avgl)
{
//TODO: Speed up..
int cnt,i;
valStruct vals[41];
double * cadata1=absarray(cdata1);
double * cadata2=absarray(cdata2);
double * cadata3=absarray(cdata3);
avgl[0]=cdata1[0];
avgl[1]=cdata2[0];
avgl[2]=cdata3[0];
// determines the 40th largest absolute value. For each color channel
valqueue vq;
cnt=0;
for (i=0;i<16384;i++) {
if (cnt==40) {
vals[40].d=cadata1[i];
vq.push(vals[40]);
vals[40]=vq.top();
vq.pop();
}
else {
vals[cnt].d=cadata1[i];
vq.push(vals[cnt]);
cnt++;
}
}
truncq(cadata1,cdata1,vq.top().d,sig1);
while(!vq.empty()) vq.pop();
cnt=0;
for (i=0;i<16384;i++) {
if (cnt==40) {
vals[40].d=cadata2[i];
vq.push(vals[40]);
vals[40]=vq.top();
vq.pop();
}
else {
vals[cnt].d=cadata2[i];
vq.push(vals[cnt]);
cnt++;
}
}
truncq(cadata2,cdata2,vq.top().d,sig2);
while(!vq.empty()) vq.pop();
cnt=0;
for (i=0;i<16384;i++) {
if (cnt==40) {
vals[40].d=cadata3[i];
vq.push(vals[40]);
vals[40]=vq.top();
vq.pop();
}
else {
vals[cnt].d=cadata3[i];
vq.push(vals[cnt]);
cnt++;
}
}
truncq(cadata3,cdata3,vq.top().d,sig3);
Safefree(cadata1);
Safefree(cadata2);
Safefree(cadata3);
return 1;
}