mandelbrot.cc

Go to the documentation of this file.

/************************************************************************/
/*                                                                      */
/*    vspline - a set of generic tools for creation and evaluation      */
/*              of uniform b-splines                                    */
/*                                                                      */
/*            Copyright 2017 - 2023 by Kay F. Jahnke                    */
/*                                                                      */
/*    Permission is hereby granted, free of charge, to any person       */
/*    obtaining a copy of this software and associated documentation    */
/*    files (the "Software"), to deal in the Software without           */
/*    restriction, including without limitation the rights to use,      */
/*    copy, modify, merge, publish, distribute, sublicense, and/or      */
/*    sell copies of the Software, and to permit persons to whom the    */
/*    Software is furnished to do so, subject to the following          */
/*    conditions:                                                       */
/*                                                                      */
/*    The above copyright notice and this permission notice shall be    */
/*    included in all copies or substantial portions of the             */
/*    Software.                                                         */
/*                                                                      */
/*    THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND    */
/*    EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES   */
/*    OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND          */
/*    NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT       */
/*    HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,      */
/*    WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING      */
/*    FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR     */
/*    OTHER DEALINGS IN THE SOFTWARE.                                   */
/*                                                                      */
/************************************************************************/
 
/// \file mandelbrot.cc
///
/// \brief  calculate an image of a section of the mandelbrot set
///
/// to demonstrate that vspline's transform routines don't have to use
/// b-splines at all, here's a simple example creating a functor
/// to perform the iteration leading to the mandelbrot set, together
/// with a vspline::domain adapting the coordinates and another
/// functor to do the 'colorization'. The three functors are chained
/// and fed to vspline::transform() to yield the image.
///
/// compile with:
///
/// ./examples.sh mandelbrot.cc
///
/// ./hwy_mandelbrot_clang++ -2 -1 1 1
///
/// ( -2 , -1 ) being lower left and ( 1 , 1 ) upper right
///
/// the result will be written to 'mandelbrot.tif'
 
#include <iostream>
#include <vspline/vspline.h>
 
#include <vigra/stdimage.hxx>
#include <vigra/imageinfo.hxx>
#include <vigra/impex.hxx>
 
// we want a colour image
typedef vigra::RGBValue<double,0,1,2> pixel_type; 
 
// coordinate_type has a 2D coordinate
typedef vigra::TinyVector < double , 2 > coordinate_type ;
 
typedef typename vspline::vector_traits < coordinate_type > :: type crd_v ;
typedef typename vspline::vector_traits < double > :: type compo_v ;
const int VSZ = vspline::vector_traits < coordinate_type > :: size ;
typedef typename vspline::vector_traits < int , VSZ > :: type int_v ;
 
// target_type is a 2D array of pixels  
typedef vigra::MultiArray < 2 , pixel_type > target_type ;
 
struct mandelbrot_functor
: public vspline::unary_functor < coordinate_type , int >
{
  typedef typename vspline::unary_functor
                   < coordinate_type , int > base_type ;
  
  using base_type::vsize ;
  
  const int max_iterations = 255 ;
  const double threshold = 1000.0 ;
  
  // the single-value evaluation recognizably uses the well-known
  // iteration formula
 
  void eval ( const coordinate_type & c , int & m ) const
  {
    std::complex < double > cc ( c[0] , c[1] ) ;
    std::complex < double > z ( 0.0 , 0.0 ) ;
    
    for ( m = 0 ; m < max_iterations ; m++ )
    {
      z = z * z + cc ;
      if ( std::abs ( z ) > threshold )
        break ;
    }
  }
  
  // the vector code is a bit more involved, since the vectorized type
  // for a std::complex<double> is a vigra::TinyVector of two SIMD types,
  // and we implement the complex maths 'manually':
 
  void eval ( const crd_v & c , int_v & m ) const
  {
    // state of the iteration
    crd_v z { 0.0f , 0.0f } ;
    // iteration count
    m = 0.0f ;
    
    for ( int i = 0 ; i < max_iterations ; i++ )
    {
      // z = z * z ; using complex maths
      compo_v rr = z[0] * z[0] ;
      compo_v ii = z[1] * z[1] ;
      compo_v ri = z[0] * z[1] ;
      z[0] = rr - ii ;
      z[1] = 2.0f * ri ;
      
      // create a mask for those values which haven't exceeded the
      // theshold yet
      
      rr += ii ;
      auto mm = ( rr < threshold * threshold ) ;
      
      // if the mask is empty, all values have exceeded the threshold
      // and we end the iteration
 
      if ( none_of ( mm ) )
        break ;
      
      // now we add 'c', the coordinate
      z += c ;
      
      // and increase the iteration count for all values which haven't
      // exceeded the threshold
      m ( mm ) = ( m + 1 ) ;
    }
  }
  
// #endif
 
} ;
 
struct colorize
: public vspline::unary_functor < int , pixel_type , VSZ >
{
  // to 'colorize' we produce black-and-white from the incoming
  // value's LSB
  
  template < class IN , class OUT >
  void eval ( const IN & c , OUT & result ) const
  {
    result = 255 * ( c & 1 ) ;
    // for a bit more colour, try
    // result = OUT( 255 * ( c & 1 ) , 127 * ( c & 2 ) , 63 * ( c & 4 ) ) ;
  } ;
} ;
 
int main ( int argc , char * argv[] )
{
  // get the extent of the section to show
  if ( argc < 5 )
  {
    std::cerr << "please pass x0, y0, x1 and y1 on the command line" << std::endl ;
    std::cerr << "typical invocation: mandelbrot -2 -1 1 1" << std::endl ;
    exit ( -1 ) ;
  }
  
  double x0 = atof ( argv[1] ) ;
  double y0 = atof ( argv[2] ) ;
  double x1 = atof ( argv[3] ) ;
  double y1 = atof ( argv[4] ) ;
  
  // this is where the result should go - we produce a 4K image:
 
  target_type target ( vigra::Shape2 ( 3840 , 2160 ) ) ;
 
  // the domain maps the image coordinates to the coordinates of the
  // section we want to display. The mapped coordinates, now in the range
  // of ((x0,y0), (x1,y1)), are fed to the functor calculating the result
  // of the iteration, and it's results are fed to a 'colorize' object
  // which translates the iteration depth values to a pixel value.
 
  auto f =   vspline::domain < coordinate_type , VSZ >
                      ( coordinate_type ( 0 , 0 ) ,
                        coordinate_type ( 3839 , 2159 ) ,
                        coordinate_type ( x0 , y0 ) ,
                        coordinate_type ( x1 , y1 ) )
           + mandelbrot_functor()
           + colorize() ;
  
  // the combined functor is passed to transform(), which uses it for
  // every coordinate pair in 'target' and deposits the result at the
  // corresponding location.
 
  vspline::transform ( f , target ) ;
  
  // store the result with vigra impex
  
  vigra::ImageExportInfo imageInfo ( "mandelbrot.tif" );
  
  std::cout << "storing the target image as 'mandelbrot.tif'" << std::endl ;
  
  vigra::exportImage ( target ,
                       imageInfo
                       .setPixelType("UINT8")
                       .setCompression("100")
                       .setForcedRangeMapping ( 0 , 255 , 0 , 255 ) ) ;
}