psychlops/extension/prototype/figure/psychlops_figure_gabor.cpp

   1 /*
   2  *  psychlops_figure_gabor.cpp
   3  *  Psychlops Standard Library (Universal)
   4  *
   5  *  Last Modified 2006/01/19 by Kenchi HOSOKAWA
   6  *  (C) 2006 Kenchi HOSOKAWA, Kazushi MARUYA and Takao SATO
   7  */
   8
   9
  10 #include <math.h>
  11 #include <vector>
  12 #include <algorithm>
  13
  14 #include "../../../platform/gl/psychlops_g_GL_h.h"
  15
  16 #include "../../../psychlops_core.h"
  17 #include "psychlops_figure_gabor.h"
  18
  19
  20 namespace Psychlops {
  21
  22
  23         double bellcurve(const double x, const double sigma) {
  24                 return exp( -(x*x) / (2.0*sigma*sigma) );
  25         }
  26
  27
  28 namespace StimulusTemplate {
  29
  30         Gabor::~Gabor() {
  31         }
  32         Gabor& Gabor::setContrast(double cont) {
  33                 contrast = cont;
  34                 return *this;
  35         }
  36         Gabor& Gabor::setWavelength(Length length) {
  37                 wavelength = length;
  38                 return *this;
  39         }
  40         Gabor& Gabor::setOrientation(Angle orient) {
  41                 orientation = orient;
  42                 return *this;
  43         }
  44         Gabor& Gabor::setPhase(Angle phs) {
  45                 phase = phs;
  46                 return *this;
  47         }
  48         const Point Gabor::getDatum() const {
  49                 Point p(left_+width_/2.0, top_+height_/2.0);
  50                 return p;
  51         }
  52
  53 }
  54
  55 namespace Figures {
  56
  57         int substructImages(Image &result, const Image &s1, const Image &s2, const double factor)
  58         {
  59                 if(s1.getHeight()!=s2.getHeight()) throw new Exception("substructImages");
  60                 if(s1.getWidth()!=s2.getWidth()) throw new Exception("substructImages");
  61                 result.release();
  62                 result.set(s1.getWidth(), s1.getHeight());
  63                 Color c1, c2, r1;
  64                 double difference = 0, sum = 0;
  65                 bool error = false, serror = false;
  66                 for(int y=0; y<s1.getHeight(); y++) {
  67                         for(int x=0; x<s1.getWidth(); x++) {
  68                                 c1 = s1.getPix(x,y);
  69                                 c2 = s2.getPix(x,y);
  70                                 difference = c1.getR()-c2.getR();
  71                                 sum += difference!=0.0 ? 1 : 0;
  72                                 error = difference>=2.0/255.0;
  73                                 serror =  serror || error;
  74                                 result.pix_raw(x,y,Color(
  75                                         error ? 1 : 0,
  76                                         difference>0 ? factor*difference : 0,
  77                                         difference<0 ? factor*difference : 0,
  78                                         //c2.getR()-c1.getR()> 0 ? factor*(c2.getR()-c1.getR()) : 0,//.5 + factor * (c1.getR() - c2.getR()),
  79                                         //c1.getR()-c2.getR()> 0 ? factor*(c1.getR()-c2.getR()) : 0,//.5 + factor * (c1.getG() - c2.getG()),
  80                                         0,//.5 + factor * (c1.getB() - c2.getB()),
  81                                         1
  82                                 ));
  83                                 //if(c1.getR()-c2.getR()>0.01) std::cout << c1.getR()-c2.getR() << " " << c1.getR() << " " <<  c2.getR() << std::endl;
  84                         }
  85                 }
  86                 return serror ? -sum : sum;
  87         }
  88
  89         void drawImageCoordinateTuner(Image &img)
  90         {
  91                 const int width=256, height=256;
  92                 img.release();
  93                 img.set(width, height);
  94                 for(int y=0; y<height; y++) {
  95                         for(int x=0; x<width; x++) {
  96                                 img.pix_raw(x,y,
  97                                   Color(abs(x-127.5)/255.0, 0.0, 0.0)
  98                                 );
  99                         }
 100                 }
 101         }
 102         void drawGrating(Image &img, int width , int height, double wavelength, double contrast, double orientation, double phase)
 103         {
 104                 img.release();
 105                 img.set(width, height);
 106                 double xp, yp, r, freq = (1.0/wavelength)*2*PI;
 107                 for(int y=0; y<height; y++) {
 108                         yp = y-height/2.0;
 109                         for(int x=0; x<width; x++) {
 110                                 xp = x-width/2.0;
 111                                 r = sqrt(xp*xp+yp*yp);
 112                                 img.pix_raw(x,y,
 113                                   Color(
 114                                         0.5+0.5 * contrast* cos(freq* (sin(orientation)*xp-cos(orientation)*yp) + phase)
 115                                   )
 116                                 );
 117                         }
 118                 }
 119         }
 120         void drawGaussian(Image &img, double sigma, double factor)
 121         {
 122                 img.release();
 123                 img.set(sigma*8, sigma*8);
 124                 const int width = img.getWidth(), height = img.getHeight();
 125                 double xp, yp, r, r2;
 126                 for(int y=0; y<height; y++) {
 127                         yp = y-height/2.0;
 128                         for(int x=0; x<width; x++) {
 129                                 xp = x-width/2.0;
 130                                 r = sqrt(xp*xp+yp*yp);
 131                                 r2 = -(r*r) / (2.0*sigma*sigma);
 132                                 img.pix_raw(x,y,Color(factor*(exp(r2))));
 133                         }
 134                 }
 135         }
 136         void drawGabor(Image &img, double sigma, double wavelength, double contrast, double orientation, double phase)
 137         {
 138                 img.release();
 139                 img.set(sigma*8, sigma*8);
 140                 const int width = img.getWidth(), height = img.getHeight();
 141                 double xp, yp, r, freq = (1.0/wavelength)*2*PI;
 142                 for(int y=0; y<height; y++) {
 143                         yp = y-height/2.0;
 144                         for(int x=0; x<width; x++) {
 145                                 xp = x-width/2.0;
 146                                 r = sqrt(xp*xp+yp*yp);
 147                                 img.pix_raw(x,y,
 148                                   Color(0.996078 * (
 149                                         0.5+0.5 * contrast* cos(freq* (sin(orientation)*xp-cos(orientation)*yp) + phase)
 150                                         * exp( -(r*r) / (2.0*sigma*sigma) ))
 151                                   ));
 152                         }
 153                 }
 154         }
 155
 156
 157
 158
 159         GaborBase::GaborBase()
 160         {
 161                 set(0,0);
 162                 contrast = 1.0;
 163                 wavelength = 10;
 164                 orientation = 0;
 165                 phase = 0;
 166         }
 167
 168         GaborBase& GaborBase::setSigma(double sigma)
 169         {
 170                 return setSigma(sigma*pixel);
 171         }
 172         GaborBase& GaborBase::setSigma(Length sigma)
 173         {
 174                 Rectangle::set(sigma*8.0, sigma*8.0);
 175                 return *this;
 176         }
 177         GaborBase& GaborBase::setWave(double wavelen, double cont, double orient, double phs)
 178         {
 179                 return setWave(wavelen*pixel, cont, orient*degree, phs*degree);
 180         }
 181         GaborBase& GaborBase::setWave(Length wavelen, double cont, Angle orient, Angle phs)
 182         {
 183                 contrast = cont;
 184                 wavelength = wavelen;
 185                 orientation = orient;
 186                 phase = phs;
 187                 return *this;
 188         }
 189
 190
 191         const bool ImageGabor::Key::operator <(Key rhs) const
 192         {
 193                 for(int i=0; i<5; i++)
 194                 {
 195                         if(data[i]!=rhs.data[i]) return data[i]<rhs.data[i];
 196                 }
 197                 return false;
 198         }
 199         const bool ImageGabor::Key::operator >(Key rhs) const
 200         {
 201                 for(int i=0; i<5; i++)
 202                 {
 203                         if(data[i]!=rhs.data[i]) return data[i]>rhs.data[i];
 204                 }
 205                 return false;
 206         }
 207
 208
 209         ImageGabor::ImageGabor() : GaborBase()
 210         {
 211         }
 212
 213         ImageGabor& ImageGabor::cache(DrawableWithCache &target)
 214         {
 215                 Image* tmp = new Image();
 216                 Figures::drawGabor(*tmp, this->getWidth()/8.0, 1.0/wavelength, contrast, orientation, phase);
 217                 tmp->cache(target);
 218                 Key key;
 219                 key.data[0] = this->getWidth();
 220                 key.data[1] = wavelength;
 221                 key.data[2] = contrast;
 222                 key.data[3] = orientation;
 223                 key.data[4] = phase;
 224                 cache_list[key] = tmp;
 225                 return *this;
 226         }
 227
 228         ImageGabor& ImageGabor::draw(Drawable &target)
 229         {
 230                 Key key;
 231                 key.data[0] = this->getWidth();
 232                 key.data[1] = wavelength;
 233                 key.data[2] = contrast;
 234                 key.data[3] = orientation;
 235                 key.data[4] = phase;
 236                 if(cache_list.count(key)!=0) {
 237                         cache_list[key]->centering(this->getCenter()).draw();
 238                 } else {
 239                         Figures::drawGabor(normal_, this->getWidth()/8.0, 1.0/wavelength, contrast, orientation, phase);
 240                         normal_.centering(this->getCenter()).draw();
 241                 }
 242                 return *this;
 243         }
 244
 245         void ImageGabor::to(Image &dest, Canvas &media)
 246         {
 247                 Figures::drawGabor(dest, this->getWidth()/8.0, 1.0/wavelength, contrast, orientation, phase);
 248
 249         }
 250
 251
 252
 253
 254 }
 255
 256         QuickGabor::Instance::Instance() : instantiated(false) {
 257                 referenced_count_ = 0;
 258         }
 259         QuickGabor::Instance::~Instance() {
 260         }
 261         void QuickGabor::Instance::release() {
 262                 referenced_count_--;
 263                 if(referenced_count_<=0) {
 264                         if(instantiated) {
 265                                 delete envelope_;
 266                                 delete [] carrier_;
 267                         }
 268                         instantiated = false;
 269                 }
 270         }
 271         void QuickGabor::Instance::set(double sigma) {
 272                 sigma_ = sigma;
 273                 size_ = (int)(sigma_*6.0);
 274                 setEnvelope();
 275                 setCarrier();
 276                 instantiated = true;
 277         }
 278         void QuickGabor::Instance::setEnvelope() {
 279                 Color col;
 280                 int particle=0;
 281                 envelope_ = new Image;
 282                 envelope_->set(size_, size_, Image::RGBA, Image::FLOAT);
 283                 int halfsize = size_/2;
 284
 285                 if(halfsize*2 != size_) particle = 1;
 286                 for(int y=-halfsize; y<halfsize+particle; y++) {
 287                         for(int x=-halfsize; x<halfsize+particle; x++) {
 288                                 //col.set( 0.5, 0.5, 0.5, 1-bellcurve(sqrt(x*x+y*y), sigma_) );
 289                                 envelope_->pix(x+halfsize, y+halfsize, Color::gray);
 290                                 envelope_->alpha(x+halfsize, y+halfsize, 1-bellcurve(::sqrt((double)(x*x+y*y)), sigma_));
 291                         }
 292                 }
 293                 envelope_->cache();
 294         }
 295         void QuickGabor::Instance::setCarrier() {
 296                 double height, halfwidth;
 297                 carrier_ = new Rectangle[size_];
 298                 for(int x=1; x<size_-1; x++) {
 299                         halfwidth = size_/2-1;
 300                         height = 2*sqrt(halfwidth*halfwidth-(x-halfwidth)*(x-halfwidth));
 301                         carrier_[x].set(1.1,height).shift(x-halfwidth-0.5,-height/2);
 302                 }
 303         }
 304         QuickGabor::Instance& QuickGabor::Instance::draw(int x, int y, Length wavelength, Angle orientation, Angle phase, double contrast, Drawable &target) {
 305                 Color col;
 306                 int limit = size_/2 -1;
 307                 glPushMatrix();
 308 //              glEnable(GL_LINE_SMOOTH);
 309 //              glHint(GL_LINE_SMOOTH_HINT, GL_NICEST);
 310                 glTranslated(x+limit,y+limit,0);
 311                 glRotated(orientation,0,0,1);
 312                 for(int i=1; i<size_-3; i++) {
 313                         col.set( contrast *
 314                                                 sin( ((double)i/(double)wavelength)*2*PI + phase.at_radian() )
 315                                         /2.0+0.5);
 316                         carrier_[i].draw(col, target);
 317                 }
 318                 glPopMatrix();
 319                 envelope_->draw(x,y,target);
 320                 return *this;
 321         }
 322
 323
 324
 325         std::vector<QuickGabor::Instance *> QuickGabor::instance_;
 326         int QuickGabor::setInstance(double sigma) {
 327                 Instance *ins;
 328                 int index = 0;
 329                 bool found =false;
 330                 if(instance_.empty()) {
 331 //                      instance_ = std::vector<Instance>();
 332                 } else {
 333                         for(std::vector<Instance *>::iterator iter=instance_.begin(); iter!=instance_.end(); iter++, index++) {
 334                                 if((*iter)->sigma_ == sigma) {
 335                                         (*iter)->referenced_count_++;
 336                                         found = true;
 337                                         break;
 338                                 }
 339                         }
 340                 }
 341                 if(found==true) {
 342                         return index;
 343                 } else {
 344                         ins = new Instance;
 345                         ins->set(sigma);
 346                         instance_.push_back( ins );
 347                         instance_.back()->referenced_count_++;
 348                         return index;
 349                 }
 350                 return -1;
 351         }
 352
 353         void QuickGabor::release() {
 354                 if(sigma_index_>=0 && sigma_index_<instance_.size()) instance_[sigma_index_]->release();
 355         }
 356
 357
 358
 359         QuickGabor::QuickGabor() {
 360         }
 361         QuickGabor::QuickGabor(double length, double sigma, double cont, double orient, double phs) {
 362                 set(length, sigma, cont, orient, phs);
 363         }
 364         QuickGabor::QuickGabor(Length length, Length sigma, double cont, Angle orient, Angle phs) {
 365                 set(length, sigma, cont, orient, phs);
 366         }
 367         QuickGabor::~QuickGabor() {
 368                 release();
 369         }
 370         QuickGabor & QuickGabor::set(double length, double sigma, double cont, double orient, double phs) {
 371                 return set(length*pixel, sigma*pixel, cont, orient*degree, phs*degree);
 372         }
 373         QuickGabor & QuickGabor::set(Length length, Length sigma, double cont, Angle orient, Angle phs) {
 374                 sigma_ = sigma;
 375                 contrast = cont;
 376                 wavelength = length;
 377                 orientation = orient;
 378                 phase = phs;
 379                 left_ = 0.0;
 380                 top_ = 0.0;
 381
 382                 width_ = (int)(sigma_*6.0);
 383                 height_ = width_;
 384
 385                 sigma_index_ = setInstance(sigma_);
 386
 387                 return *this;
 388         }
 389         QuickGabor & QuickGabor::setDatum(const Point &po) {
 390                 return centering(po);
 391         }
 392         QuickGabor & QuickGabor::shift(const double x, const double y, const double z) {
 393                 left_ += x;
 394                 top_ += y;
 395                 return *this;
 396         }
 397         QuickGabor & QuickGabor::centering(const double x, const double y, const double z) {
 398                 left_ = x - width_/2.0;
 399                 top_ = y - height_/2.0;
 400                 return *this;
 401         }
 402         QuickGabor & QuickGabor::centering(const Point &po) {
 403                 return centering(po.x, po.y);
 404         }
 405         QuickGabor & QuickGabor::centering(const Drawable &target) {
 406                 return centering(Display::getCenter());
 407         }
 408
 409         QuickGabor & QuickGabor::draw(Drawable &target) {
 410                 instance_[sigma_index_]->draw((int)left_, (int)top_, wavelength, orientation, phase, contrast, target);
 411                 return *this;
 412         }
 413         QuickGabor & QuickGabor::draw(int x, int y, Drawable &target) {
 414                 instance_[sigma_index_]->draw(x, y, wavelength, orientation, phase, contrast, target);
 415                 return *this;
 416         }
 417
 418 }
 419