/* CLIP_A */
/*




		CLIP 1.04



		C++ Library for Image Processing
		Version 1.04
		May 2005

		Filename: "picture.h"

		John A Robinson
		http://www.intuac.com/userport/john/clip104.html




  This program is free software; you can redistribute it and/or modify
  it under the terms of the GNU General Public Licence as published by
  the Free Software Foundation; either version 2 of the Licence, 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 Licence for more details.
  
  Portions of this software are copyright as follows:
	Extracts from the CImg class library
		Copyright David Tschumperlé
		GNU General Public Licence
	The Small JPEG decoder library
		Copyright Rich Geldreich, 1994 - 2000
		GNU Lesser General Public Licence
	Inverse DCT functions from The Independent JPEG Group
		Copyright Tom Lane, 1991 - 1998
		Bespoke licence
	Extracts from OpenIllusionist
		Copyright Justen Hyde and Dan Parnham, 2004.
		GNU General Public Licence
	CameraFramework class
		Copyright Enrico Costanza, 2003.
		GNU General Public Licence
	Adaptive Prediction Trees encoding and decoding
		Copyright John Robinson 1995 - 2004
		Bespoke licence
	Extracts from the Numerics library
		Copyright Brent Worden, 1998 - 1999.
		Bespoke licence (unrestricted use)
	Extracts from the Bioinformatics Template Library
		Copyright Birkbeck College, 1997 - 1998.
		GNU General Public Licence
	All other code
		Copyright John Robinson 1994 - 2005
		GNU General Public Licence

  The bespoke licences referred to above all allow free use and
  adaptation of code, but some (e.g. the GNU Lesser Public Licence 
  and the APT licence) allow use in non-free packages, whereas
  the GNU General Public Licence requires that any adapted
  versions are also free and available in source form. This
  program is released under the General Public Licence because
  it uses other code that was released under that licence.

  Contents:

	Major section		Starting tag (search for this to
				find the start of the section)
	--------------------------------------------------------------
	Commentary header			CLIP_A
	Header for doxygen			CLIP_B
	Start of code				CLIP_C
	Declarations of the CLIP types		CLIP_D
		This documents the API and
		has doxygen-formatted 
		commentary. It is equivalent
		to the old CLIP header file
		when the rest of the package
		was a library of separate
		source files.
	Image display 				CLIP_E
	Camera input 				CLIP_F
	JPEG decoding 				CLIP_G
	APT encoding/decoding 			CLIP_H
	Generic image file input/output		CLIP_I
	Font and text picture constructors	CLIP_J
	CLIP core functions			CLIP_K
	CLIP matrix functions			CLIP_L

	

*/
/* CLIP_B */
/// \file picture.h
///
/// CLIP is a C++ class library for manipulating pictures with integer
/// pel values. It has a simple programming model, efficient support for 
/// common image processing operations, Matlab-like constructs for 
/// subimages, real-time capture and display. It works on MS Windows,
/// Linux and Macintosh platforms.
///
/// CLIP is distributed as a single header file that incorporates all
/// functionality. I recommend that you use precompilation of the
/// header if your compiler supports it.

/* CLIP_C */
#ifndef _CLIP_1_0_1
#define _CLIP_1_0_1

#include <cstdio>
#include <cstdlib>
#include <cstdarg>
#include <cfloat>
#include <cmath>
#include <cstring>
#include <ctime>
#include <iostream>
#include <fstream>
#include <cctype>
#include <setjmp.h>
#include <assert.h>
#include <limits.h>


// Put the following inside the clip scope if you don't like std being
// included in global.
using namespace std;



/* CLIP_D */

// Begin effective header file: All CLIP types are declared in this section.
// CLIP classes are in the clip namespace. But if you look at the end of
// this file, you'll see the directive: using namespace clip; . This might
// be evil, but it means that old CLIP programs work without change.
// Start of CLIP scope
namespace clip {

// Have to declare things from clip's internals that are used in the
// type declarations.
namespace clip_internal {
	struct CImgDisplay;
	template <typename PEL> class double_dereferenced_picture_of_;
	}
// And the classes everyone is going to be friends with:
	template <typename PEL>  class picture_of_;
	template <typename PEL>  class colour_picture_of_;


/// Represents a range of integers with constant spacing, e.g. 2, 5, 8, .., 32

/// iranges are used chiefly to index picture_of_ints and colour_picture_of_ints 
/// to define subimages.
class irange {
	int start, end, spacing;
	int num_in_irange;
public:

/// irange r(0,1,7); means r represents the range of ints 
/// 0,1,2,3,4,5,6,7.  Similarly, irange s(3,7,31); means
/// s represents 3,7,11,15,19,23,27,31.
    irange(const int& first, const int& second, const int& last);

/// irange r(5); means that the range is just the single integer 5.
    irange(const int& solevalue);

/// Conventional copy constructor
    irange(const irange& other);

/// Conventional assignment.
    irange& operator= (const irange& other);

/// Instantiates range containing just the integer 0
    irange();

/// Destructor
    ~irange();

/// Returns number of integers in the range
    int length();

/// Example: If we declared irange r(0,2,10); then the numbers in 
/// the range are 0,2,4,6,8,10 and r[3] equals 6.
    int operator[](const int& i) const;
 
/// Test for ranges identical
    int operator== 	(const irange& other) const;

/// Test for ranges different
    int operator!= (const irange& other) const;

/// Example: irange r(0,2,10); r+=3; Then r is 3,5,7,9,11,13.
    irange& operator+= (const int& i);

/// Other inc and dec operators are analogous to +=
    irange& operator-= (const int& i);
    irange& operator++ ();
    irange& operator-- ();
    irange operator++ (int);
    irange operator-- (int);

	// The following functions are used internally in CLIP to avoid
	// excessive friendship (which VC6 doesn't seem to like).
	// They are probably not useful in applications.
	int length() const {
		return num_in_irange;
		}
	int istart() const {
		return start;
		}
	int iend() const {
		return end;
		}
	int ispacing() const {
		return spacing;
		}		
    
    };


//
//  Declarations for global functions that use iranges
//
irange operator+(const irange& o, const int& i);
irange operator+(const int& i, const irange& o);
irange operator-(const irange& o, const int& i);
irange operator-(const int& i, const irange& o);

//--------------------------------------
// picture_of_int and ancilliary classes
//--------------------------------------


// The pic_rep class is used internally by CLIP 
namespace clip_internal {
template <typename PEL> class pic_rep{
	PEL **buf;
	int border;
	int rows, cols;
	int numrefs;
	pic_rep(const int& r, const int& c, const int& bord) {
	int i;		// Counter
	rows = r;
	cols = c;
        border = bord;
        buf = new PEL *[rows+2*border];
	buf[0] = new PEL[(rows+2*border)*(cols+2*border)];
	buf[0] += border;
        for (i = 1; i < rows+2*border; i++)
                {
		buf[i] = buf[0] + i*(cols+2*border);
                }
        buf += border;
        // Now set the border to zero
        for (i = -border; i < 0; i++)
                memset(buf[i]-border,0,(cols+2*border)*sizeof(PEL));
        for (; i < rows; i++) {
                memset(buf[i]-border,0,border*sizeof(PEL));
                memset(buf[i]+cols,0,border*sizeof(PEL));
                }
        for (; i < rows+border; i++)
                memset(buf[i]-border,0,(cols+2*border)*sizeof(PEL));
        }
template <typename OTHERPEL>
// See below for case when OTHERPEL is same type as PEL
	pic_rep(const pic_rep<OTHERPEL>& other)
	{
	// cerr << "pic rep copy constr\n";
	rows = other.rows;
	cols = other.cols;
        border = other.border;
        buf = new PEL *[rows+2*border];
	buf[0] = new PEL[(rows+2*border)*(cols+2*border)];
	buf[0] += border;
	int i;
        for (i = 1; i < rows+2*border; i++)
                {
		buf[i] = buf[0] + i*(cols+2*border);
                }
        buf += border;
	// Now copy data over
	for (i = -border; i < rows+border; i++)
	for (int j = -border; j < cols+border; j++)
		buf[i][j] = (PEL)other.buf[i][j];
	}
// template <>
	pic_rep(const pic_rep<PEL>& other)
	{
	// cerr << "pic rep copy constr\n";
	rows = other.rows;
	cols = other.cols;
        border = other.border;
        buf = new PEL *[rows+2*border];
	buf[0] = new PEL[(rows+2*border)*(cols+2*border)];
	buf[0] += border;
        for (int i = 1; i < rows+2*border; i++)
                {
		buf[i] = buf[0] + i*(cols+2*border);
                }
        buf += border;
	// Now copy data over
	memcpy(&buf[-border][-border],&other.buf[-border][-border],
			(rows+2*border)*(cols+2*border)*sizeof(PEL));
	}

	void dc_pad()
	{
	int i, j;
	for (i = 0; i < rows; i++) {
		for (j = -border; j < 0; j++)
			buf[i][j] = buf[i][0];
		for (j = cols; j < cols+border; j++)
			buf[i][j] = buf[i][cols-1];
		}
	for (i = -border; i < 0; i++)
		memcpy(buf[i]-border,buf[0]-border,(cols+2*border)*sizeof(PEL));
	for (i = rows; i < rows+border; i++)
		memcpy(buf[i]-border,buf[rows-1]-border,
			(cols+2*border)*sizeof(PEL));
	}
	void zero_border()
	{
	int i;
	for (i = 0; i < rows; i++) {
		memset(buf[i]-border,0,border*sizeof(PEL));
		memset(buf[i]+cols,0,border*sizeof(PEL));
		}
	for (i = -border; i < 0; i++)
		memset(buf[i]-border,0,(cols+2*border)*sizeof(PEL));
	for (i = rows; i < rows+border; i++)
		memset(buf[i]-border,0,
			(cols+2*border)*sizeof(PEL));
	}
	~pic_rep() {
		delete [] (buf[-border]-border);
        buf -= border;
        delete [] buf;
	}
// Have to deal with friends differently because of MSC VC++ problem
#if ( defined(_MSC_VER)) 
	// Instead have to do
	friend class picture_of_<unsigned char>;
	friend class picture_of_<char>;
	friend class picture_of_<unsigned short>;
	friend class picture_of_<short>;
	friend class picture_of_<unsigned int>;
	friend class picture_of_<int>;
	friend class picture_of_<unsigned long>;
	friend class picture_of_<long>;
	friend class picture_of_<float>;
	friend class picture_of_<double>;
	friend class picture_of_<long double>;
	friend class double_dereferenced_picture_of_<unsigned char>;
	friend class double_dereferenced_picture_of_<char>;
	friend class double_dereferenced_picture_of_<unsigned short>;
	friend class double_dereferenced_picture_of_<short>;
	friend class double_dereferenced_picture_of_<unsigned int>;
	friend class double_dereferenced_picture_of_<int>;
	friend class double_dereferenced_picture_of_<unsigned long>;
	friend class double_dereferenced_picture_of_<long>;
	friend class double_dereferenced_picture_of_<float>;
	friend class double_dereferenced_picture_of_<double>;
	friend class double_dereferenced_picture_of_<long double>;
	friend class colour_picture_of_<unsigned char>;
	friend class colour_picture_of_<char>;
	friend class colour_picture_of_<unsigned short>;
	friend class colour_picture_of_<short>;
	friend class colour_picture_of_<unsigned int>;
	friend class colour_picture_of_<int>;
	friend class colour_picture_of_<unsigned long>;
	friend class colour_picture_of_<long>;
	friend class colour_picture_of_<float>;
	friend class colour_picture_of_<double>;
	friend class colour_picture_of_<long double>;
	friend class pic_rep<unsigned char>;
	friend class pic_rep<char>;
	friend class pic_rep<unsigned short>;
	friend class pic_rep<short>;
	friend class pic_rep<unsigned int>;
	friend class pic_rep<int>;
	friend class pic_rep<unsigned long>;
	friend class pic_rep<long>;
	friend class pic_rep<float>;
	friend class pic_rep<double>;
	friend class pic_rep<long double>;
#else
	template <class APEL> friend class picture_of_;
	template <class APEL> friend class double_dereferenced_picture_of_;
	template <class APEL> friend class colour_picture_of_;
	template <class APEL> friend class pic_rep;
#endif
	};
}

const int DEFAULT_BORDER = 32; 		// When a picture is instantiated 
					// without an explicit border size, 
					// this is the size used. (See the 
					// picture_of_int constructors below).

/// The main class in CLIP. A picture_of_ is a one-plane (greylevel) image.

/// picture_of_s are the basic type in CLIP: colour_picture_of_s are made out
/// of them and for greylevel processing, they are all you need.
template <typename PEL> class picture_of_ {
// Private part of picture_of_int
    // Following data items are used for sequence processing:
    bool seq;
    int numfiles;
    int currfile;
    bool usingcam;
    int firstcycle;
    char **av;
    picture_of_<PEL> *backup;

// Have to deal with friends differently because of MSC VC++ problem
#if (defined(_MSC_VER)) 
	friend class picture_of_<unsigned char>;
	friend class picture_of_<char>;
	friend class picture_of_<unsigned short>;
	friend class picture_of_<short>;
	friend class picture_of_<unsigned int>;
	friend class picture_of_<int>;
	friend class picture_of_<unsigned long>;
	friend class picture_of_<long>;
	friend class picture_of_<float>;
	friend class picture_of_<double>;
	friend class picture_of_<long double>;
	friend class clip_internal::double_dereferenced_picture_of_<unsigned char>;
	friend class clip_internal::double_dereferenced_picture_of_<char>;
	friend class clip_internal::double_dereferenced_picture_of_<unsigned short>;
	friend class clip_internal::double_dereferenced_picture_of_<short>;
	friend class clip_internal::double_dereferenced_picture_of_<unsigned int>;
	friend class clip_internal::double_dereferenced_picture_of_<int>;
	friend class clip_internal::double_dereferenced_picture_of_<unsigned long>;
	friend class clip_internal::double_dereferenced_picture_of_<long>;
	friend class clip_internal::double_dereferenced_picture_of_<float>;
	friend class clip_internal::double_dereferenced_picture_of_<double>;
	friend class clip_internal::double_dereferenced_picture_of_<long double>;
	friend class colour_picture_of_<unsigned char>;
	friend class colour_picture_of_<char>;
	friend class colour_picture_of_<unsigned short>;
	friend class colour_picture_of_<short>;
	friend class colour_picture_of_<unsigned int>;
	friend class colour_picture_of_<int>;
	friend class colour_picture_of_<unsigned long>;
	friend class colour_picture_of_<long>;
	friend class colour_picture_of_<float>;
	friend class colour_picture_of_<double>;
	friend class colour_picture_of_<long double>;
#else
	template <class APEL> friend class picture_of_;
	template <class APEL> friend class clip_internal::double_dereferenced_picture_of_;
	template <class APEL> friend class colour_picture_of_;
#endif
protected:
	clip_internal::pic_rep<PEL> *pr;
	struct clip_internal::CImgDisplay *qxp;
	irange rowrange, colrange;
	int dereferenced;	// 1 = singly, 2 = doubly

// Function for averaging three channels into one (colour to monochrome)
template <typename OTHERPEL>
	picture_of_& average_assign(const picture_of_<OTHERPEL>& rb, 
	    const picture_of_<OTHERPEL>& rc, const picture_of_<OTHERPEL>& rd) {
        if ((rowrange.length() != rb.rowrange.length())||
            (colrange.length() != rb.colrange.length()))
                cerr << "warning: range["<<rowrange.length()<<"]["<<
                        colrange.length()<<"]=range["<<
                        rb.rowrange.length()<<"]["<<
                        rb.colrange.length()<<"]\n";
        int outr = rowrange.istart();
        int inr = rb.rowrange.istart();
        while((outr <= rowrange.iend())&&(inr <= rb.rowrange.iend()))
                {
                row_average_op(pr->buf[outr],rb.pr->buf[inr],
			rc.pr->buf[inr],rd.pr->buf[inr],colrange,rb.colrange);
                outr += rowrange.ispacing();
                inr += rb.rowrange.ispacing();
                }
        return *this;
	}

// Templated for dealing with all kinds of update between all kinds of
// picture_of objects
template <typename OTHERPEL>
	picture_of_& assign(const picture_of_<OTHERPEL>& rb, const char& sel) {
        if ((rowrange.length() != rb.rowrange.length())||
            (colrange.length() != rb.colrange.length()))
                cerr << "warning: range["<<rowrange.length()<<"]["<<
                        colrange.length()<<"]=range["<<
                        rb.rowrange.length()<<"]["<<
                        rb.colrange.length()<<"]\n";
        int outr = rowrange.istart();
        int inr = rb.rowrange.istart();
        while((outr <= rowrange.iend())&&(inr <= rb.rowrange.iend()))
                {
                row_update_op(sel,pr->buf[outr],rb.pr->buf[inr],
                        colrange, rb.colrange);
                outr += rowrange.ispacing();
                inr += rb.rowrange.ispacing();
                }
        return *this;
        }

// Non-templated for dealing with assignment of same kind of picture
	picture_of_& assign(const picture_of_<PEL>& rb) {
        if ((rowrange.length() != rb.rowrange.length())||
            (colrange.length() != rb.colrange.length()))
                cerr << "warning: range["<<rowrange.length()<<"]["<<
                        colrange.length()<<"]=range["<<
                        rb.rowrange.length()<<"]["<<
                        rb.colrange.length()<<"]\n";
        int outr = rowrange.istart();
        int inr = rb.rowrange.istart();
        while((outr <= rowrange.iend())&&(inr <= rb.rowrange.iend()))
                {
                row_assign_op(pr->buf[outr],rb.pr->buf[inr],
                        colrange, rb.colrange);
                outr += rowrange.ispacing();
                inr += rb.rowrange.ispacing();
                }
        return *this;
        }

	picture_of_& assign(const PEL& newval, const char& sel) {
	// Fill templine, then do other lines by copying
	PEL *temp = new PEL[colrange.length()];
	for (int j = 0; j < colrange.length(); j++)
		temp[j] = newval;
	irange temprange(0,1,colrange.length()-1);
	for (int i = rowrange.istart(); i <= rowrange.iend(); i += rowrange.ispacing())
		{
		if ((sel == '=')&&(colrange.ispacing() == 1)) // Do fast
			memcpy(&pr->buf[i][colrange.istart()],temp,
				colrange.length()*sizeof(PEL));
		else
			row_update_op(sel, pr->buf[i], temp, 
				colrange, temprange);
		}
	delete [] temp;
	return *this;
	}


	PEL sum(const picture_of_<PEL>& rb, const char& sel) const {
        if ((rowrange.length() != rb.rowrange.length())||
            (colrange.length() != rb.colrange.length()))
                cerr << "warning: picture["<<rowrange.length()<<"]["<<
                        colrange.length()<<
                        "] combined and summed with picture["<<
                        rb.rowrange.length()<<"]["<<
                        rb.colrange.length()<<"]\n";
        PEL sum = 0;
        int outr = rowrange.istart();
        int inr = rb.rowrange.istart();
        while((outr <= rowrange.iend())&&(inr <= rb.rowrange.iend()))
                {
                sum += row_sum_op(sel,pr->buf[outr],rb.pr->buf[inr],
			colrange, rb.colrange);
                outr += rowrange.ispacing();
                inr += rb.rowrange.ispacing();
                }
        return sum;
        }

	void row_mask_op(const int& thresh, const int& below,
    		const int& equal, const int& above, 
		register PEL *dest, register PEL *source,
		const irange& destrange, const irange& sourcerange){
	register int i = destrange.length();
	if (sourcerange.length() < i)
		i = sourcerange.length();
	dest += destrange.istart();
	source += sourcerange.istart();
	register int dinc = destrange.ispacing();
	register int sinc = sourcerange.ispacing();
	register PEL current;
	for (;i--;dest+=dinc,source+=sinc) {
		current = *source;
		if (current > thresh)
			*dest = above;
		else if (current == thresh)
			*dest = equal;
		else
			*dest = below;
		}
	}

	void row_limits_op(const char& sel, const register PEL *source,
		const irange& sourcerange, PEL& currlim) const {
	register int i = sourcerange.length();
	source += sourcerange.istart();
	register int sinc = sourcerange.ispacing();
	register PEL lim = currlim;
	switch (sel) {
		case 'n':	// Find minimum
			for (; i--; source+=sinc) {
				if (*source < lim) 
					lim = *source;
				}
			break;
		case 'x':	// Find maximum
			for (; i--; source+=sinc) {
				if (*source > lim) 
					lim = *source;
				}
			break;
		}
	currlim = lim;
	}

	void row_limits_op(const char& sel, register PEL *source,
		const irange& sourcerange, const PEL& newlim) {
	register int i = sourcerange.length();
	source += sourcerange.istart();
	register int sinc = sourcerange.ispacing();
	switch (sel) {
		case 'N':	// Set minimum (clip)
			for (; i--; source+=sinc) {
				if (*source < newlim) 
					*source = newlim;
				}
			break;
		case 'X':	// Set maximum (clip)
			for (; i--; source+=sinc) {
				if (*source > newlim) 
					*source = newlim;
				}
			break;
		}
	}
		

	void row_assign_op(register PEL *dest, register PEL *source, 
		const irange& destrange, const irange& sourcerange) {
	register int i = destrange.length();
	if (sourcerange.length() < i)
		i = sourcerange.length();
	dest += destrange.istart();
	source += sourcerange.istart();
	register int dinc = destrange.ispacing();
	register int sinc = sourcerange.ispacing();
	if ((sinc==1)&&(dinc==1))
		memcpy(dest,source,i*sizeof(PEL));
	else 	{
		for (;i--;dest+=dinc,source+=sinc)
			*dest = *source;
		}
	}

template <typename OTHERPEL>
	void row_average_op(
		register PEL *dest, register OTHERPEL *s1, 
		register OTHERPEL *s2, register OTHERPEL *s3, 
		const irange& destrange, const irange& sourcerange) {
	register int i = destrange.length();
	if (sourcerange.length() < i)
		i = sourcerange.length();
	dest += destrange.istart();
	s1 += sourcerange.istart();
	s2 += sourcerange.istart();
	s3 += sourcerange.istart();
	register int dinc = destrange.ispacing();
	register int sinc = sourcerange.ispacing();
		for (;i--;dest+=dinc,s1+=sinc,s2+=sinc,s3+=sinc)
			*dest = (PEL)((*s1+*s2+*s3)/3);
	}

template <typename OTHERPEL>
	void row_update_op(const char& sel, 
		register PEL *dest, register OTHERPEL *source, 
		const irange& destrange, const irange& sourcerange) {
	register int i = destrange.length();
	if (sourcerange.length() < i)
		i = sourcerange.length();
	dest += destrange.istart();
	source += sourcerange.istart();
	register int dinc = destrange.ispacing();
	register int sinc = sourcerange.ispacing();
	switch (sel) {
		case '=':
			for (;i--;dest+=dinc,source+=sinc)
				*dest = (PEL) *source;
			break;
		case '+':
			for (;i--;dest+=dinc,source+=sinc)
				*dest += (PEL) *source;
			break;
		case '-':
			for (;i--;dest+=dinc,source+=sinc)
				*dest -= (PEL) *source;
			break;
		case '*':
			for (;i--;dest+=dinc,source+=sinc)
				*dest *= (PEL) *source;
			break;
		case '/':
			for (;i--;dest+=dinc,source+=sinc) {
				if (*source == 0) {
					cout << "Attempt to divide by a pel ";
					cout << "with value 0. Dividing by ";
					cout << "1 instead." << endl;
					}
				else
					*dest /= (PEL) *source;
				}
			break;
		case '|':	// Absolute difference
			for (;i--;dest+=dinc,source+=sinc)
				{
				register PEL tempd = *dest;
				register PEL temps = (PEL) *source;
				if (tempd > temps)
					*dest = tempd - temps;
				else
					*dest = temps - tempd;
				}
			break;
		case '^':	// Squared difference
			for (;i--;dest+=dinc,source+=sinc)
				{
				register PEL temp = *dest - (PEL) *source;
				*dest = temp*temp;
				}
			break;
		case '&':	// Conditional replacement
			for (;i--;dest+=dinc,source+=sinc)
				{
				register PEL zero = 0;
				// Above line to avoid compiler warnings
				// when comparing unsigned types.
				register PEL temp = (PEL) *source;
				if (temp >= zero)
					*dest = temp;
				}
			break;
		default:
			for (;i--;dest+=dinc,source+=sinc)
				*dest = (PEL) *source;
			break;
		}
	}
	PEL row_sum_op(const char& sel, 
		register PEL *first, register PEL *second, 
		const irange& firstrange, const irange& secondrange) const {
	register int i = firstrange.length();
	if (secondrange.length() < i)
		i = secondrange.length();
	first += firstrange.istart();
	second += secondrange.istart();
	register int finc = firstrange.ispacing();
	register int sinc = secondrange.ispacing();
	PEL sum = 0;
	switch (sel) {
		case '=':
			for (;i--;first+=finc,second+=sinc)
				sum += *second;
			break;
		case '+':
			for (;i--;first+=finc,second+=sinc)
				sum += *first + *second;
			break;
		case '-':
			for (;i--;first+=finc,second+=sinc)
				sum += *first - *second;
			break;
		case '*':
			for (;i--;first+=finc,second+=sinc)
				sum += *first * *second;
			break;
		case '/':
			for (;i--;first+=finc,second+=sinc) {
				if (*second == 0) {
					cout << "Attempt to divide by a pel ";
					cout << "with value 0. Dividing by ";
					cout << "1 instead." << endl;
					sum += *first;
					}
				else
					sum += *first / *second;
				}
			break;
		case '|':	// Absolute difference
			for (;i--;first+=finc,second+=sinc)
				{
				register PEL temp = *first - *second;
				if (temp < 0)
					temp = -temp;
				sum += temp;
				}
			break;
		case '^':	// Squared difference
			for (;i--;first+=finc,second+=sinc)
				{
				register PEL temp = *first - *second;
				sum += temp*temp;
				}
			break;
		default:
			for (;i--;first+=finc,second+=sinc)
				sum += *second;
			break;
		}
	return sum;
	}

	void init_new(const int& r, const int& c, const int& bord) {
		using clip_internal::pic_rep;
		irange temprow(0,1,r-1);
		irange tempcol(0,1,c-1);
		rowrange = temprow;
		colrange = tempcol;
        	dereferenced = 0;
		qxp = 0;	// Initialized off
		pr = new pic_rep<PEL>(r, c, bord);
        }
	void fake_new(const picture_of_<PEL>& other) {
// Used only for setting up colour pictures that are pure monochrome
		const int r = other.nrows();
		const int c = other.ncols();
		irange temprow(0,1,r-1);
		irange tempcol(0,1,c-1);
		rowrange = temprow;
		colrange = tempcol;
        	dereferenced = 0;
		qxp = 0;	// Initialized off
		pr = other.pr;
        }
// Following constructor sets up dereferenced version
	picture_of_(const irange& range, const picture_of_<PEL>& source)
	// The only contructor that doesn't use init_new()
        {
		pr = source.pr;	// Note: no copying of data
		qxp = 0;	// Initialized off
        if (source.dereferenced == 0)
                {
                dereferenced = 1;
                rowrange = range;
                colrange = source.colrange;
                }
        else if (source.dereferenced == 1)
                {
                dereferenced = 2;
                rowrange = source.rowrange;
                colrange = range;
                }
        }
// And an equivalent function called from colour_picture_of_<PEL>
	void dereffrom(const irange& range, const picture_of_<PEL>& source)
	// Exact copy of the above constructor
	// Called only after null constructor in colour_picture_of_
        {
		pr = source.pr;	// Note: no copying of data
		qxp = 0;	// Initialized off
        if (source.dereferenced == 0)
                {
                dereferenced = 1;
                rowrange = range;
                colrange = source.colrange;
                }
        else if (source.dereferenced == 1)
                {
                dereferenced = 2;
                rowrange = source.rowrange;
                colrange = range;
                }
        }
//
// Following constructor and updater are only accessed by the
// friend class colour_picture_of_
//
	picture_of_<PEL>(void) {
		// Does nothing
	}
	void initialize(const int r, const int c, const int bord = DEFAULT_BORDER) {
        init_new(r,c,bord);
        for (int i = 0; i < r; i++)
                memset(pr->buf[i],0,c*sizeof(PEL));
        }
	void fake_initialize(const picture_of_<PEL>& other) {
// Called from a colour picture that is pure monochrome
        fake_new(other);
        }

	int construct_or_read_from_file(const char *pnmname,
		const int bord = DEFAULT_BORDER, const int readflag = 1);

// Now follows the part accessible by applications
public:
/// Instantiates a picture by reading it in from a file. E.g.
/// picture_of_int in("goldhill256"); sets up in and reads in goldhill256.
/// The following file formats are read: PGM, PPM, JPEG, APT and the most
/// common BMP formats.
/// If the file is in colour, it is read and converted to monochrome.
/// The bord argument specifies how many pels border to add around the
/// picture. You don't have to allow for the border when indexing the
/// picture. Having a border allows you to do, e.g.,  in[-1][-1]
/// without causing a segmentation fault. This is useful when writing
/// local operators that slide over the picture: You don't have to treat
/// the edges specially.
    picture_of_(const char *pnmname, const int bord=DEFAULT_BORDER) {
	if (construct_or_read_from_file(pnmname, bord, 0) >= 0)
		return;
	qxp = 0;
	char errmsg[180];
	sprintf(errmsg,"CLIP cannot read\n%s\nas a picture",pnmname);
	picture_of_<PEL> temp(errmsg,0,255);
	int r=temp.nrows();
	int c=temp.ncols();
        init_new(r,c,bord);
        for (int i = 0; i < r; i++)
        for (int j = 0; j < c; j++)
                pr->buf[i][j] = temp[i][j];
	}


/// Instantiates a zeroed pic of given dimensions. E.g. 
/// picture_of_int pic(240,320); creates a zeroed (black) picture	
/// of size 320 columns by 240 rows.	
    picture_of_(const int r, const int c, const int bord=DEFAULT_BORDER) {
        init_new(r,c,bord);
        for (int i = 0; i < r; i++)
                memset(pr->buf[i],0,c*sizeof(PEL));
        }
	
/// Creates a picture from a 1D array of the pels in scan order. Rarely used.
    picture_of_(const int r, const int c, const PEL *in, 
		const int bord=DEFAULT_BORDER) {	
        init_new(r,c,bord);
        for (int i = 0; i < r; i++)
                {
                memcpy(pr->buf[i],in,c*sizeof(PEL));
                in += c;
                }
        }

/// Conventional copy constructor
template <typename OTHERPEL>
    picture_of_(const picture_of_<OTHERPEL>& other) {
	using clip_internal::pic_rep;
	dereferenced = other.dereferenced;
	rowrange = other.rowrange;
	colrange = other.colrange;
	qxp = 0;	// Initialized off
	if (dereferenced)
		// Copy constructor does not copy pic_rep
		pr = (pic_rep<PEL> *) other.pr;
		// Above cast ok because never dereference to another type.
	else
		pr = new pic_rep<PEL>(*other.pr);
	}
// And the special case of constructing from same type that VC doesn't
// seem to handle as part of above constructor.
// template <>
    picture_of_(const picture_of_<PEL>& other) {
	using clip_internal::pic_rep;
	dereferenced = other.dereferenced;
	rowrange = other.rowrange;
	colrange = other.colrange;
	qxp = 0;	// Initialized off
	if (dereferenced)
		// Copy constructor does not copy pic_rep
		pr = other.pr;
		// Above cast ok because never dereference to another type.
	else
		pr = new pic_rep<PEL>(*other.pr);
	}


/// Create from a colour_picture_of_int (includes conversion to monochrome)
template <typename OTHERPEL>
    picture_of_(const colour_picture_of_<OTHERPEL>& other)
	{
	using clip_internal::pic_rep;
	qxp = 0;	// Initialized off
	int colismono = other.ismonochrome();
	if ((colismono && other.mono.dereferenced)||
	    (!colismono && other.red.dereferenced)) {
		cerr << "Sorry. You cannot construct a picture_of_ ";
		cerr << "from an\nirange-dereferenced colour_picture_of_. ";
		cerr << "Use assignment.\n";
		exit(1);
		}
	dereferenced = 0;
	if (colismono) {
		rowrange = other.mono.rowrange;
		colrange = other.mono.colrange;
		pr = new pic_rep<PEL>(*other.mono.pr);
		}
	else {
		rowrange = other.red.rowrange;
		colrange = other.red.colrange;
		int r = rowrange.length();
		int c = colrange.length();
        	int bord = other.red.pr->border;
        	init_new(r,c,bord);
		average_assign(other.green,other.red,other.blue);
		}
	}

/// Modify-source constructor that applies callback point function (*f)
/// over source to create new picture.
    picture_of_(const picture_of_<PEL>& source, PEL (*f)(PEL&, const PEL&)) {
        int r = source.rowrange.length();
        int c = source.colrange.length();
        int bord = source.pr->border;
        init_new(r, c, bord);
		PEL *p;
		int outrow = 0;
		for (int i = source.rowrange.istart(); i <= source.rowrange.iend();
	    i += source.rowrange.ispacing())
		{
		p = pr->buf[outrow++];
		for (int j = source.colrange.istart(); j <= source.colrange.iend();
		    j += source.colrange.ispacing())
			f(*p++, source.pr->buf[i][j]);
		}
	}

/// Modify-source constructor that applies callback neighbourhood function (*f)
/// over source to create new picture.
		picture_of_(const picture_of_<PEL>& source, int (*f)(PEL&, const PEL **)) {
		int i,j;	// Counters
        int r = source.rowrange.length();
        int c = source.colrange.length();
        int bord = source.pr->border;
        init_new(r, c, bord);
        // Have to reset pr row pointers so f works as expected
        int totalrows = source.pr->rows;
        int sourcecoloff = source.colrange.istart();
        for (i = 0-bord; i < totalrows+bord; i++)
                source.pr->buf[i] += sourcecoloff;
        for (j = 0; j < c; j++)
            {
            int sourcerowoff = source.rowrange.istart();
            for (i=0; i<r; i++,sourcerowoff+=source.rowrange.ispacing())
                f(pr->buf[i][j], (const PEL **) &source.pr->buf[sourcerowoff]);
		// Throw away the return val because this is a constructor
            for (i = 0-bord; i < totalrows+bord; i++)
                source.pr->buf[i] += source.colrange.ispacing();
            sourcecoloff += source.colrange.ispacing();
            }
        for (i = 0-bord; i < totalrows+bord; i++)
                source.pr->buf[i] -= sourcecoloff;
        }

/// The ONLY way to create text (e.g. labels) in CLIP is to instantiate
/// a picture_of_int with the appropriate text. This creates a picture with
/// lettering level fg on a ground level bg. size is an integer scale.
/// The result can be combined into your pictures as required.
    picture_of_(const char *text, const int fg, const int bg,	
			const int size=1);		

/// Read picture from file into object. The picture_of_int object is NOT
/// resized, so read only reads as much of the input file as will fit.
    int read(const char *name);		
					
/// Write picture to a PGM file. PGM is CLIP's "native" file format.		
    int write(const char *pgmname);	

/// Write picture in APT format. 0 <= quality <= 100. 100=lossless.
    int aptwrite(const char *aptname, const int quality = 100);

/// Write picture in BMP format. If you must.
    int bmpwrite(const char *bmpname);

/// Load from a colour_picture_of_, stacking into a colour format that can
/// be used by functions like showascolour().
    int loadascolour(const colour_picture_of_<PEL>& other);


/// Destructor
    ~picture_of_() {
	// cerr << "picture_of_ delete\n";
	if (dereferenced)
		return;
	if (qxp) 
		delete qxp;
	delete pr;
        }


/// Set values in the border of the picture by copying the values on
/// the edges of the picture outwards.
    void dc_pad() { pr->dc_pad(); }	 

/// Set values in the border of the picture to zero.
    void zero_border() { pr->zero_border(); }

/// Reflects picture at boundaries to prevent edge effects for a filter of
/// horizontal length hlength and vertical length vlength.
/// When convolving an image whose left hand side is
/// a b c d e f ...
/// with a filter with an odd number of taps, the picture should be 
/// 'reflected' in the boundary as:
/// ... f e d c b | a b c d e f ...
/// where | represents the boundary.
/// Convolving with a filter with an even number of taps, requires reflecting:
/// ... f e d c b a | a b c d e f ...
/// symmetric_extension does the appropriate reflection based on given lengths.
    void symmetric_extension(int hlength, int vlength) {
	int heven = ((hlength>>1)<<1 == hlength); // = 1 if even num of taps
	int veven = ((vlength>>1)<<1 == vlength); // = 1 if even num of taps
	for (int i = 0; i < nrows(); i++) {
		for (int j = -hlength/2; j < 0; j++) {
			pr->buf[i][j] = pr->buf[i][-j-heven];
			pr->buf[i][ncols()-1-j] = pr->buf[i][ncols()-1+j+heven];
			}
		}
	for (int j = 0; j < ncols(); j++) {
		for (int i = -vlength/2; i < 0; i++) {
			pr->buf[i][j] = pr->buf[-i-veven][j];
			pr->buf[nrows()-1-i][j] = pr->buf[nrows()-1+i+veven][j];
			}
		}
	}	

					 
/// Return number of rows in picture				
    int nrows()  const			
	{return rowrange.length();}

/// Return number of columns in picture
    int ncols()  const			
	{return colrange.length();}

/// Return width of border
    int bordersize() const
	{return pr->border; }

/// Return maximum value in picture and write its position into 
/// reti (row) and retj (column).
    PEL max(int& reti, int& retj) const {
	PEL retval = pr->buf[0][0];
	for (int i = rowrange.istart(); i <= rowrange.iend(); 
	    i += rowrange.ispacing())
		{
		for (int j = colrange.istart(); j <= colrange.iend(); 
		    j += colrange.ispacing())
			if (pr->buf[i][j] > retval)
				{
				retval = pr->buf[i][j];
				reti = i;
				retj = j;
				}
		}
	return retval;
	}

					
/// Return minimum value in picture and write its position into 
/// reti (row) and retj (column).
    PEL min(int& reti, int& retj) const {
	PEL retval = pr->buf[0][0];
	for (int i = rowrange.istart(); i <= rowrange.iend(); 
	    i += rowrange.ispacing())
		{
		for (int j = colrange.istart(); j <= colrange.iend(); 
		    j += colrange.ispacing())
			if (pr->buf[i][j] < retval)
				{
				retval = pr->buf[i][j];
				reti = i;
				retj = j;
				}
		}
	return retval;
	}


/// Return maximum value in picture
    PEL max() const {
	PEL retval = pr->buf[0][0];
	for (int i = rowrange.istart(); i <= rowrange.iend(); 
	    i += rowrange.ispacing())
		row_limits_op('x',pr->buf[i],colrange,retval);
	return retval;
	}

/// Return minimum value in picture
    PEL min() const {
	PEL retval = pr->buf[0][0];
	for (int i = rowrange.istart(); i <= rowrange.iend(); 
	    i += rowrange.ispacing())
		row_limits_op('n',pr->buf[i],colrange,retval);
	return retval;
	}

/// Clip picture at given maximum value
    void max(const PEL clipval) {
	for (int i = rowrange.istart(); i <= rowrange.iend(); 
	    i += rowrange.ispacing())
		row_limits_op('X',pr->buf[i],colrange,clipval);
	}

/// Clip picture at given minimum value
    void min(const PEL clipval) {
	for (int i = rowrange.istart(); i <= rowrange.iend(); 
	    i += rowrange.ispacing())
		row_limits_op('N',pr->buf[i],colrange,clipval);
	}

/// Return total of all pel values
    PEL total() const {
	PEL retval = 0;
	for (int i = rowrange.istart(); i <= rowrange.iend(); 
	    i += rowrange.ispacing())
		{
		for (int j = colrange.istart(); j <= colrange.iend(); 
		    j += colrange.ispacing())
			retval += pr->buf[i][j];
		}
	return retval;
	}

/// \name Sequence processing facilities
/// Sequence processing facilities for "standard" CLIP programs.
/// Many CLIP programs consist of a processing loop that cycles over
/// input pictures from a sequence of named files or from a camera.

//@{

/// Sequence processing constructor that interprets main()'s parameters. 
/// If the command line had no arguments, tries to attach a camera. 
/// If that fails, generates a default picture. 
/// If the command line had arguments, interprets these as picture
/// files and sets up for reading these in sequence.
/// See also member function next()
picture_of_<PEL>(int argc, char *argv[], 
		const int def_nrows=240, const int def_ncols=320, 
		const int bord=DEFAULT_BORDER) {
	seq = 1;	// Flag to say sequence processing
	numfiles = argc - 1;
	currfile = 1;
	usingcam = 0;
	firstcycle = 1;
	av = argv;
	qxp = 0;
	backup = 0;
	if (!numfiles) {
		// If there are no arguments, try to open a camera and, if
		// that fails, create a default image
		initialize(def_nrows, def_ncols, bord);
		if (attach_camera() < 0) {
			hramp();
			clip::picture_of_<PEL> text("Default image",0,128,3);
			paste_centred_on(text, def_ncols/2,def_nrows/2);
			}
		else {
			request_frame();
			usingcam = 1;
			}
		return;
		}
	// Otherwise, open the first file
	if (construct_or_read_from_file(argv[1], bord, 0) >= 0) {
		seq = 1;
		if (numfiles == 1)
			// Make and save a backup for fast reread
			backup = new picture_of_(*this);
		return;
		}
	qxp = 0;
	char errmsg[180];
	sprintf(errmsg,"CLIP cannot read\n%s\nas a picture",argv[1]);
	clip::picture_of_<PEL> temp(errmsg,0,255);
	int r=temp.nrows();
	int c=temp.ncols();
	initialize(r, c, bord);
	*this = temp;
	if (numfiles == 1)
		// Make and save a backup for fast reread
		backup = new picture_of_(*this);
	}

/// Loads the next picture in the sequence into this picture that was
/// constructed with the sequence processing constructor. Returns 1 if
/// the picture has been seen before (if cycling back to the start of
/// the argument list) and 0 if it has not. Returns -1 if the picture
/// was not constructed using the sequence processing constructor.
int next()
	{
	if (!seq)	// Not set up to do sequence processing
		{
		cerr << "next() called for a picture not constructed for ";
		cerr << "sequence processing." << endl;
		return -1;
		}
	if (usingcam) {
		while(!frame_delivered())
			;
		request_frame();
		return 1;	// Not seen this picture before
		}
	if (!numfiles){		// Using default picture
		// Remake it, just in case it was overwritten
		hramp();
		clip::picture_of_<PEL> text("Default image",0,128,3);
		paste_centred_on(text,ncols()/2,nrows()/2);
		return 0;
		}
	if (numfiles == 1){	// We made a backup of this one to save
				// rereading every cycle
		*this = *backup;
		return 0;
		}
	if (currfile > numfiles) {
		currfile = 1;	// Back to start
		firstcycle = 0;
		}
	read(av[currfile++]);
	return firstcycle;
	}

//@}


/// Reranges or normalizes picture values to fall in given range
/// If scale_even_if_inside_range is false (default), then if all
/// values are current between the give newlow and newhigh, they
/// are unaltered. Otherwise, values are scaled so that the new
/// max() will be newhigh and the new min() will be newmin.
    int rerange(const PEL newlow, const PEL newhigh, 
    			const bool scale_even_if_inside_range = 0) {
    	if (newlow == newhigh) {
		*this = newhigh;
		return 0;
		}
    	PEL mx = max();
	PEL mn = min();
	if ((mx == newhigh) && (mn == newlow))
		return 0;
	if (!scale_even_if_inside_range && (mx <= newhigh) && (mn >= newlow))
		return 0;
	if (mx == mn) {
		if (mx > 0)	*this = newhigh;
		else		*this = newlow;
		return 0;
		}
	PEL newrange = newhigh - newlow;
	PEL oldrange = mx - mn;
	for (int i = 0; i < nrows(); i++) {
		PEL *pel = pr->buf[i];
		for (int j = 0; j < ncols(); j++) {
			*pel = newlow + (newrange*(*pel-mn))/oldrange;
			pel++;
			}
		}
	return 1;
	}

/// \name Picture display functions
/// Use these functions to show or reshow the current picture on the screen

//@{

/// Display the picture on the screen in its own window, with
/// the comment string as title. If normalize is true values are scaled
/// to [0,255] for display. If false, it is assumed they already are in this range.
    int show(const char *comment, const bool normalize=1);

/// Display picture without a specific title.
    int show(const bool n = 1)
	{ return show("CLIP image",n); }

/// Display without decorations at Nx,Ny. If parent==0, the coords are
/// screen coordinates, and the window receives events (e.g. pointx())
/// as normal. If parent is a displayed picture, the coords are relative to
/// that picture and a child window is displayed that does not receive events.
    int showat(const int Nx, const int Ny, const bool normalize,
			const picture_of_ *parent); 

/// Display without decorations, with colour_picture_of_int as parent
    int showat(const int Nx, const int Ny, const bool normalize,
	const colour_picture_of_<PEL> *parent); 	

/// Display without decorations at screen coords Nx, Ny
	int showat(const int Nx, const int Ny, const bool normalize) 
		{ return showat(Nx,Ny,normalize,(const picture_of_ *)0); }

/// Display without decorations at screen coords Nx, Ny
	int showat(const int Nx, const int Ny) 
		{ return showat(Nx,Ny,1,(const picture_of_ *)0); }

/// Redisplay the picture (presumably because it has changed) and change
/// title to the comment string.
    int reshow(const char *comment,	const bool normalize=1);

/// Redisplay the picture without altering the title.
    int reshow(const bool n = 1);	

/// Change title of the displayed pic.
    int rename(const char *comment);
	
/// Stop showing the picture
    int unshow();

/// Display a picture_of_int as if each integer pel were a
/// colour vector of four unsigned chars in the order alpha (most
/// significant) - red - green - blue.  Use comment string as window title.
	int showascolour(const char *comment);	

/// As above, without window title
    int showascolour()
		{ return showascolour("CLIP image"); }

/// Redisplay a picture_of_int as if each integer pel were a colour
/// vector of four unsigned chars in the order alpha-red-green-blue.
/// Change title to comment string.
    int reshowascolour (const char *comment); 		

/// As above, without change in window title
    int reshowascolour();

/// Display a picture_of_ with showascolour() interpretation (see
/// above), but without decorations at Nx,Ny. If parent==0, the coords are
/// screen coordinates, and the window receives events (e.g. pointx())
/// as normal. If parent is a displayed picture, the coords are relative to
/// that picture and a child window is displayed that does not receive events.
/// As showascolour() but wihout decorations.    
	int showascolourat(const int Nx, const int Ny,
		const picture_of_ *parent);

/// As above with colour_picture_of_ parent
    int showascolourat(const int Nx, const int Ny,
		const colour_picture_of_<PEL> *parent); 

/// As above with no parent, i.e. at screen coordinates Nx, Ny
	int showascolourat(const int Nx, const int Ny) 
		{ return showascolourat(Nx,Ny,(const picture_of_ *)0); }


/// Returns non-zero if the user has left-clicked in the close box of the
/// window where the picture is being displayed.
    int closerequested();		

/// Return current pointer (mouse) x coord relative to left of picture window    
	int pointx();

/// Return current pointer (mouse) y coord relative to top of picture window    
    int pointy();

/// Return non-zero if the mouse's left button is pressed within the picture's window.
    int pointleftbutton();
	
/// Returns non-zero if mid button down
	int pointmiddlebutton();
	
/// Returns non-zero if right button down
    int pointrightbutton();	
	
/// Returns non-zero if the shift key is held down while a mouse button is down.
    int pointshiftbutton();		

/// Returns keycode of any key held down    
	int key();

/// Clears keystroke so that it will not be reprocessed.
    int clearkey();	

/// Inspect a shown picture. Until the right button is clicked inside
/// the picture, the program pauses and allows the user to
/// left-click within the picture and inspect the values at that point.	
    int inspect() {	
	if (!qxp) return -1;
	char infostring[64];
	int prevx = -1;
	int prevy = -1;
	while (pointrightbutton());	// From any previous rightbutton press
	while(1) {
		if (pointrightbutton()) {
			// No longer do a rename back to what was before
			return 1;
			}
		if (!pointleftbutton())
			continue;
		int x = pointx();
		int y = pointy();
		if ((x == prevx)&&(y == prevy))
			continue;
		prevx = x; prevy = y;
		if ((x < 0)||(x >= ncols())||(y < 0)||(y >= nrows())) {
			// Beyond picture limits
			rename("pointer outside pic");
			continue;		// Back to start of loop
			}
		// Else pointer is within picture limits
		float pelvalue = (float) pr->buf[y][x];
		sprintf(infostring,"value at (%d,%d) = %f", x,y,pelvalue);
		rename(infostring);	// Update header
		}
	}
	
//@}

/// \name Dereferencing operators
/// The following operator[]s are used to dereference (i.e. index) pictures.
/// They must always be appended to a picture_of_int in pairs, with the
/// first [] specifying a row value or range and the second [] specifying
/// a column value or range. E.g.: picture_of_int a(r,c); a[y][x] = 255;
/// So when using []s, picture_of_int looks like a two dimensional array.
/// The difference is it can be indexed by iranges and doubles as well as
/// by ints.

//@{


/// When dereferenced (i.e. indexed) with irange operators, a subimage or
/// block of the picture is returned as a picture. E.g. irange rg(0,1,7);
/// picture_of_ pic(256,256); pic[rg][rg] = 255; sets the top left 8x8 square 
/// in pic to white.
    picture_of_<PEL>  operator[]	(const irange& selector) const {
        picture_of_<PEL> m(selector,*this);
        return m; }


/// When dereferenced by an integer, just gives entry into pel buffer.
/// No bounds checking. Same as rowptr(selector);
    PEL * operator[]  (const int& selector)
		{ return pr->buf[selector]; }
    PEL * operator[]  (const int& selector) const
		{ return pr->buf[selector]; }

/// When dereferenced (i.e. indexed) with a double, a PEL is returned which
/// is the value bilinearly interpolated at the given double coords. E.g.
/// val = a[20.5][31.435]; sets val to a weighted average of a[20][31],
/// a[20][32], a[21][31] and a[21][32]. You cannot use double indexing on
/// the left hand side of an expression because the pel at the given location
/// doesn't really exist. I.e. saying a[20.5][31.435] = val; is illegal.
/// (This works via the helper class double_dereferenced_picture_of_int
/// which users do not access directly.)
    clip_internal::double_dereferenced_picture_of_<PEL> operator[](const double selector) {	
        clip_internal::double_dereferenced_picture_of_<PEL> m(selector, this);
        return m; }

/// We often want to paste one picture on another, centred at a particular
/// location. Can do this with irange derefencings, but this way is easier.
    int paste_centred_on(const picture_of_<PEL>& from, const int x, const int y) {
	irange r(0,1,from.nrows()-1);
	irange c(0,1,from.ncols()-1);
	r -= from.nrows()/2;
	c -= from.ncols()/2;
	r += y;
	c += x;
	operator[](r).operator[](c).assign(from,'=');
	return 0;
	}

/// As above, for colour_picture_of_ pasting
    int paste_centred_on(const colour_picture_of_<PEL>& from, const int x, const int y) {
	irange r(0,1,from.nrows()-1);
	irange c(0,1,from.ncols()-1);
	r -= from.nrows()/2;
	c -= from.ncols()/2;
	r += y;
	c += x;
	operator[](r).operator[](c).assign(from,'=');
	return 0;
	}


//@}

/// The map operators provide a fast way of mapping a bilinearly interpolated
/// value from one picture into another. transparency is a number from 0 (overwrite) 
/// to 255 (preserve almost all of what was there originally). So with two
/// picture_of_int's source and dest, dest.map(50,76,source,5.3,99.99,0);
/// is exactly equivalent to dest[76][50] = source[99.99][5.3]; but much faster.
    int map(const int myx,const int myy, const picture_of_& from,	
		const double fromx, const double fromy,	
    	const int transparency);

/// This map operator is even faster than the double version above. This
/// time fromx and fromy are both 256 times the actual values desired, so
/// dest.map(50,76,source,64,384,0); is the same as 
/// dest.map(50,76,source,0.25,1.5,0); but faster. Rarely used.
     int map(const int myx,const int myy,
		const picture_of_& from, const int fromx,const int fromy,
    	const int transparency);

/// \name point and neighbourhood operators via callbacks
//@{

/// Apply callback point function *f() over picture, modifying itself.
    PEL point(PEL (*f)(PEL&))
	{
	PEL accum = 0;
	for (int r = rowrange.istart(); r <= rowrange.iend(); r += rowrange.ispacing())
	for (int c = colrange.istart(); c <= colrange.iend(); c += colrange.ispacing())
		accum += f(pr->buf[r][c]);
	return accum;
	}

/// Apply callback point function *f() over source picture rb, putting
/// result into *this.
    PEL point(const picture_of_<PEL>& rb, PEL (*f)(PEL&, const PEL&))	
	{
        if ((rowrange.length() != rb.rowrange.length())||
            (colrange.length() != rb.colrange.length()))
                cerr << "warning: picture["<<rowrange.length()<<"]["<<
                        colrange.length()<<
                        "] in a point operation with picture["<<
                        rb.rowrange.length()<<"]["<<
                        rb.colrange.length()<<"]\n";
	PEL accum = 0;
	int outr, inr;
	int outc, inc;
	for (outr = rowrange.istart(), inr = rb.rowrange.istart(); 
	    ((outr <= rowrange.iend()) && (inr <= rb.rowrange.iend())); 
	    outr += rowrange.ispacing(), inr += rb.rowrange.ispacing())
	for (outc = colrange.istart(), inc = rb.colrange.istart();
	    ((outc <= colrange.iend()) && (inc <= rb.colrange.iend()));
	    outc += colrange.ispacing(), inc += rb.colrange.ispacing())
		accum += f(pr->buf[outr][outc], rb.pr->buf[inr][inc]);
	return accum;
	}

/// Apply callback neigh function *f() over source picture rb, putting
/// result into this.
    PEL neigh(picture_of_<PEL>& rb, PEL(*f)(PEL&, const PEL **))
	{
        if ((rowrange.length() != rb.rowrange.length())||
            (colrange.length() != rb.colrange.length()))
                cerr << "warning: picture["<<rowrange.length()<<"]["<<
                        colrange.length()<<
                        "] in a neighbourhood operation with picture["<<
                        rb.rowrange.length()<<"]["<<
                        rb.colrange.length()<<"]\n";
	PEL accum = 0;
	int outr, inr;
	int outc, inc;
	int i;
	// Have to reorganise rb pointers so that ops in f are always
	// relative to p[0][0]
	int totalrows = rb.pr->rows;
	int rbcoloff = rb.colrange.istart();
	for (i = 0; i < totalrows; i++)
		rb.pr->buf[i] += rbcoloff;
	for (outc = colrange.istart(), inc = rb.colrange.istart();
	    ((outc <= colrange.iend()) && (inc <= rb.colrange.iend()));
	    outc += colrange.ispacing(), inc += rb.colrange.ispacing())
		{
		for (outr = rowrange.istart(), inr = rb.rowrange.istart(); 
		    ((outr <= rowrange.iend()) && (inr <= rb.rowrange.iend())); 
		    outr += rowrange.ispacing(), inr += rb.rowrange.ispacing())
			accum += f(pr->buf[outr][outc],
				(const int **) &rb.pr->buf[inr]);
		for (i = 0; i < totalrows; i++)
			rb.pr->buf[i] += rb.colrange.ispacing();
		rbcoloff += rb.colrange.ispacing();
		}
	for (i = 0; i < totalrows; i++)
		rb.pr->buf[i] -= rbcoloff;
	
	return accum;
	}
//@}


/// Set pels in this picture, according to their values in source picture
/// rb: if source pel is below thresh, set dest pel to value below,
/// if source pel is equal to thresh, set dest pel to value equal,
/// if source pel is above thresh, set dest pel to int value above.
    picture_of_& mask		
		(const picture_of_& rb, const PEL& thresh,const PEL& below,
		const PEL& equal,  const PEL& above) {
        if ((rowrange.length() != rb.rowrange.length())||
            (colrange.length() != rb.colrange.length()))
                cerr << "warning: range["<<rowrange.length()<<"]["<<
                        colrange.length()<<"]=range["<<
                        rb.rowrange.length()<<"]["<<
                        rb.colrange.length()<<"]\n";
        int outr = rowrange.istart();
        int inr = rb.rowrange.istart();
        while((outr <= rowrange.iend())&&(inr <= rb.rowrange.iend()))
                {
                row_mask_op(thresh, below, equal, above,
			pr->buf[outr],rb.pr->buf[inr],
                        colrange, rb.colrange);
                outr += rowrange.ispacing();
                inr += rb.rowrange.ispacing();
                }
        return *this;
        }


/// Return 1 if and only if every pel in the picture equals testval.
	int operator==(PEL testval) {
	for (int i = rowrange.istart(); i <= rowrange.iend(); i += rowrange.ispacing())
		{
		PEL *p = pr->buf[i];
		p += colrange.istart();
		int j = colrange.length();
		for (; j--; p += colrange.ispacing())
			if (*p != testval)
				return 0;
		}
	return 1;
	}

					
/// \name Assignment operators
//@{
/// Set picture_of_ equal to colour_picture_of_. Does conversion to monochrome
template <typename OTHERPEL>
    picture_of_<PEL>& operator= (const colour_picture_of_<OTHERPEL>& rb) {
	// Convert a colour_picture_of_ to monochrome
	if (rb.ismonochrome()) {
		assign(rb.mono,'=');
		return *this;
		}
	average_assign(rb.green,rb.red,rb.blue);
	return *this;
	}
	
/// Set this picture equal to source picture with different data type. 
/// This and all assignments
/// work with irange-referenced pictures.
template <typename OTHERPEL>
    picture_of_& operator= (const picture_of_<OTHERPEL>& rb){	 
    	return assign(rb, '='); } 	

/// Set this picture equal to source picture with same pel data type. 
// template <>
    picture_of_& operator= (const picture_of_<PEL>& rb){	 
    	return assign(rb); } 	
	
/// Set every pel to newval. E.g.irange rg(0,8,256-8);
/// picture_of_ a(256,256); a[rg][rg] = 255;
/// sets every eighth pel horizontally and vertically to white.
    picture_of_& operator=(PEL newval){		
    	return assign(newval, '='); }


/// Add source picture to this picture. I.e. each pel value in rb is
/// added to the corresponding pel's value in this.
template <typename OTHERPEL>
    picture_of_& operator+=	(const picture_of_<OTHERPEL>& rb) {
    	return assign(rb, '+'); }

/// Add newval to every pel
    picture_of_& operator+= (PEL newval){
    	return assign(newval, '+'); }

/// Subtract source picture from this picture. I.e. each pel value
/// in rb is subtracted from the corresponding pel's value in this.
template <typename OTHERPEL>
    picture_of_& operator-=	(const picture_of_<OTHERPEL>& rb) {
    	return assign(rb, '-'); }

/// Subtract newval from every pel.
    picture_of_& operator-=(PEL newval){
    	return assign(newval, '-'); }

/// Multiply source picture with this picture. I.e. each pel value
/// in rb is multiplied with the corresponding pel's value in this.
template <typename OTHERPEL>
    picture_of_& operator*=(const picture_of_<OTHERPEL>& rb) {
    	return assign(rb, '*'); }

/// Multiply every pel by newval
    picture_of_& operator*=(PEL newval){
    	return assign(newval, '*'); }

/// If a pel in rb is 0, divides the corresponding pel in this
/// by 1, otherwise divides it by rb pel value.
template <typename OTHERPEL>
    picture_of_& operator/=(const picture_of_<OTHERPEL>& rb) {
    	return assign(rb, '/'); }

/// Divide every pel by newval (unless newval == 0)
    picture_of_& operator/=(PEL newval){
    	return assign(newval, '/'); }

/// The |= operator does a pelwise absolute difference, i.e.
/// destpel = abs(sourcepel-destpel) not a pelwise OR.
template <typename OTHERPEL>
    picture_of_& operator|=(const picture_of_<OTHERPEL>& rb) { 
    	return assign(rb, '|'); }	
 
/// Replace every pel with its absolute difference from newval
    picture_of_& operator|=(PEL newval){
    	return assign(newval, '|'); }

/// The ^= operator does a pelwise squared difference, i.e.
/// destpel = sourcepel-destpel; destpel *= destpel;
/// not a pelwise XOR.
template <typename OTHERPEL>
    picture_of_& operator^=(const picture_of_<OTHERPEL>& rb) {	
    	return assign(rb, '^'); }	
					
/// Replace every pel with its squared difference from newval
    picture_of_& operator^= (PEL newval){
    	return assign(newval, '^'); }

/// The &= operator does a pelwise conditional replacement, i.e.
/// if (sourcepel < 0) destpel=destpel; else destpel = sourcepel;
template <typename OTHERPEL>
    picture_of_& operator&=	(const picture_of_<OTHERPEL>& rb) {	
    	return assign(rb, '&'); }	
 
// Following two operators were included in earlier versions of
// clip for completeness, but they really have no value!
    picture_of_& operator&= (PEL newval){
    	return assign(newval, '&'); }

//@}

/// \name Operators combining pictures to give PEL results
//@{

/// When a picture is added to a picture the result is a PEL which is the
/// sum of all pelwise additions.
    PEL operator+ (const picture_of_& rb) {	
    	return sum(rb, '+'); }

/// Return the sum of all pelwise subtractions.
    PEL operator-(const picture_of_& rb) {
    	return sum(rb, '-'); }

/// Return the sum of all pelwise multiplications. Can convolve a picture with
/// an operator by making the operator into a picture, then multiply operator
/// and a moving block of the picture together.
    PEL operator*(const picture_of_& rb) {
    	return sum(rb, '*'); }

/// Return the sum of all pelwise divisions. Uses divisor of 1 for any zero
/// pel in rb.
    PEL operator/ (const picture_of_& rb) {
    	return sum(rb, '/'); }

/// Return the sum of all pelwise absolute differences. Can be used for block matching.
    PEL operator| (const picture_of_& rb) {
    	return sum(rb, '|'); }


/// Return the sum of all pelwise squared differences. Can be used for block matching.
    PEL operator^ (const picture_of_& rb) {
    	return sum(rb, '^'); }


//@}

/// Returns the PEL * buffer pointer used internally by picture_of_int. THIS
/// FUNCTION IS DANGEROUS. It providesno border management. Rarely used.
    PEL **bufptr() { return pr->buf; }
	
/// Returns the int buffer pointer used internally by picture_of_int for
/// the row i. Use this function if you need direct access to the pel values
/// to implement a fast row operation. 
     PEL *rowptr(int i) 			
	{ return pr->buf[i]; }
	 
/// const version of the above for read only cases.
   PEL *rowptr(int i) const
	{ return pr->buf[i]; }

/// \name Functions for interfacing to a local camera
//@{

/// If a camera exists, attach it to this picture. Returns negative value 
/// if the camera doesn't exist, could not be initialized, is already 
/// attached to another picture or the current picture is too small to 
/// have a camera attached.
	int attach_camera();

/// Initiate the request of a frame from the camera attached to this 
/// picture. Returns immediately; you must check frame_delivered() for 
/// arrival of the frame.
	int request_frame();

/// This version of request_frame() only works with Philips cameras on 
/// Linux. Will be extended to other cameras and platforms when possible.
	int request_frame(int shutterspeed, int gain);

/// Immediately returns 0 if a requested frame is not yet available. If 
/// the frame is available does a buffer switch and returns 1. Because of 
/// the buffer switch within the function, immediately before the call the 
/// picture has its old contents; immediately after it has its the frame.
	int frame_delivered();

/// Call to detach the camera from the picture. You may need to do this 
/// in order to attach the camera to a different picture.
	int detach_camera();
//@}

/// \name Miscellaneous matrix manipulation functions
//@{

/// Pour the pels of the source picture into this picture which has the
/// same number of pels but a different aspect ratio
int reaspect(const picture_of_<PEL>& source) {
	int snr = source.nrows();
	int snc = source.ncols();
	int dnr = nrows();
	int dnc = ncols();
	if ((snr == dnr)&&(snc == dnc)) {
		*this = source;
		return 0;
		}
	if (snr*snc != dnr*dnc) {
		cout << "Cannot reformat a " << snr << " x " << snc;
		cout << " image into a " << dnr << " x " << dnc;
		cout << " image.\nDifferent number of pels\n";
		return -1;
		}
	int cnt(snr*snc);
	int si(source.rowrange.istart());
        int sj(source.colrange.istart());
        int di(rowrange.istart());
        int dj(colrange.istart());
	int sisp(source.rowrange.ispacing());
        int sjsp(source.colrange.ispacing());
        int disp(rowrange.ispacing());
        int djsp(colrange.ispacing());
	int siend(source.rowrange.iend());
        int sjend(source.colrange.iend());
        int diend(rowrange.iend());
        int djend(colrange.iend());
	while(cnt--) {
		if (sj > sjend) {
			si += sisp;
			sj = source.colrange.istart();
			}
		if (dj > djend) {
			di += disp;
			dj = colrange.istart();
			}
		pr->buf[di][dj] = source[si][sj];
		dj += djsp;
		sj += sjsp;
		}
	/*
	int si(0), sj(0), di(0), dj(0);
	while(cnt--) {
		if (sj == snc) {
			si++;
			sj = 0;
			}
		if (dj == dnc) {
			di++;
			dj = 0;
			}
		pr->buf[di][dj++] = source[si][sj++];
		}
	*/
	return 0;
	}

/// Matrix multiply into this picture.
/// The appropriate picture dimensions have to be correct. To just do a
/// linear transform of a picture into another picture, see mx_transform
/// below.
int mx_mul(picture_of_<PEL>& A, picture_of_<PEL>& B) {
	int i, j, k;
	const int l(B.ncols());
	const int n(B.nrows());
	const int m(A.nrows());
	if (l != ncols()) {
		cout << "Cannot do C=A*B matrix multiply: Ccols (";
		cout << ncols() << ") != Bcols (" << l << ")\n";
		return -1;
		}
	if (n != A.ncols()) {
		cout << "Cannot do C=A*B matrix multiply: Acols (";
		cout << A.ncols() << ") != Brows (" << n << ")\n";
		return -1;
		}
	if (m != nrows()) {
		cout << "Cannot do C=A*B matrix multiply: Crows (";
		cout << nrows() << ") != Arows (" << m << ")\n";
		return -1;
		}
	for (i = 0; i < m; i++) {
		for (j = 0; j < l; j++) {
			PEL sum = 0;
			for (k = 0; k < n; k++)
				sum += A[i][k]*B[k][j];
			pr->buf[i][j] = sum;
			}
		}
	return 0;
	}
/// Apply picture transform as a linear transform to source and put the 
/// result in this.  Essentially a matrix multiplication, but source and 
/// this are treated as vectors. This is a very liberal function, allowing
/// sizes of vectors that are smaller than the transform.
int mx_transform(picture_of_<PEL>& transform, picture_of_<PEL>& source,
	const int use_transpose = 0) {
	int i, k;
	int snr = source.nrows();
	int snc = source.ncols();
	int sn = snr*snc;
	int dnr = nrows();
	int dnc = ncols();
	int dn = dnr*dnc;
	const int l(1);		// Effectively a column vector
	int n(transform.ncols());
	int m(transform.nrows());
	if (use_transpose) { i = n; n = m; m = i; }
	if (n < sn) {
		cout << "Cannot do linear transform:  mx cols < pels in source\n";
		return -1;
		}
	if (m < dn) {
		cout << "Cannot do linear transform:  mx rows < pels in dest\n";
		return -1;
		}
	// Note that it doesn't matter if the matrix is larger than the
	// number of pels in the source or dest.
	n = sn;
	m = dn;
	int si(0), sj(0), di(0), dj(0);
	// cerr << "mx_transform: " << nrows() << "x" << ncols() << " = ";
	// cerr << transform.nrows() << "x" << transform.ncols() << " * ";
	// cerr << source.nrows() << "x" << source.ncols() << endl;
	for (i = 0; i < m; i++) {
		si = 0;
		sj = 0;
		PEL sum = 0;
		if (use_transpose) {
		    for (k = 0; k < n; k++) {
			sum += transform[k][i]*source[si][sj];
			if (sj == snc-1) {
				si++;
				sj = 0;
				}
			else
				sj++;
			}
		    }
		else {
		    for (k = 0; k < n; k++) {
			sum += transform[i][k]*source[si][sj];
			if (sj == snc-1) {
				si++;
				sj = 0;
				}
			else
				sj++;
			}
		    }
		pr->buf[di][dj] = sum;
		if (dj == dnc-1) {
			di++;
			dj = 0;
			}
		else
			dj++;
		}
	return 0;
	}
	
/// Make this picture into the transpose of the source picture
int mx_transpose(picture_of_<PEL>& A) {
	int i, j;
	int l = A.nrows();
	int m = A.ncols();
	if (l != ncols()) {
		cout << "Cannot do matrix transpose: Bcols != Arows\n";
		return -1;
		}
	if (m != nrows()) {
		cout << "Cannot do matrix transpose: Brows != Acols\n";
		return -1;
		}
	for (i = 0; i < l; i++)
	for (j = 0; j < m; j++)
		pr->buf[j][i] = A[i][j];
	return 0;
	}

/// Make this picture into the matrix inverse of the source picture
int mx_inverse(picture_of_<PEL>& A);

/// Make the first row of this picture into the eigenvalues of the symmetric
/// matrix represented as source picture A.
int mx_eigenvalues(picture_of_<PEL>& A);

/// Make the first row of this picture into the eigenvalues of the symmetric
/// matrix represented as source picture A and write the corresponding 
/// eigenvectors to evecs
int mx_eigensystem(picture_of_<PEL>& A, picture_of_<PEL>& evecs);

/// Given a eigenvalue for the symmetric matrix represented as source 
/// source picture A, write the corresponding
/// eigenvector into the first row of this picture
int mx_eigenvector(picture_of_<PEL>& A, PEL evalue);

/// Solve linear equations A.x = b without inverting A. this is x; both
/// this and b should be one-row pictures.
int mx_solve(picture_of_<PEL>& transform, picture_of_<PEL>& b);

/// Return the determinant of this picture, treated as a matrix
PEL mx_determinant();

//@}

/// \name Miscellaneous graphics functions
//@{

/// Fill the picture with a horizontal ramp
	int hramp() {	
	// Fill picture with a horizontal ramp
	for (int i = 0; i < nrows(); i++)
	for (int j = 0; j < ncols(); j++)
		pr->buf[i][j] = j;
	return 0;
	}

/// Fill the picture with a vertical ramp
	int vramp() {	
	for (int i = 0; i < nrows(); i++)
	for (int j = 0; j < ncols(); j++)
		pr->buf[i][j] = i;
	return 0;
	}

/// Draw line of shade value from (x0,y0) to (x1,y1)
	int draw_line(int x0, int y0, int x1, int y1, PEL value) {
	int dx = x1-x0;
	int dy = y1-y0;
	if ((x0 >= 0)&&(x1 >= 0)&&(y0 >= 0)&&(y1 >= 0)&&(x0 < ncols())&&
		(x1 < ncols())&&(y0 < nrows())&&(y1 < nrows())) {
		// No clipping needed
	    if( abs(dx)>abs(dy) ){
		if( x0<x1 ){
			for(int x=x0;x<x1;x++){
				int y = y0 + dy * (x-x0) / dx;
				pr->buf[y][x] = value;
			}
		}else{
			for(int x=x0;x>x1;x--){
				int y = y0 + dy * (x-x0) / dx;
				pr->buf[y][x] = value;
			}
		}
	    }else{
		if( y0<y1 ){
			for(int y=y0;y<y1;y++){
				int x = x0 + dx * (y-y0) / dy;
				pr->buf[y][x] = value;
			}
		}else{
			for(int y=y0;y>y1;y--){
				int x = x0 + dx * (y-y0) / dy;
				pr->buf[y][x] = value;
				}
			}
		}
	    }
	else {
	    if( abs(dx)>abs(dy) ){
		if( x0<x1 ){
			for(int x=x0;x<x1;x++){
				if (x < 0 || x >= ncols())
					continue;
				int y = y0 + dy * (x-x0) / dx;
				if (y < 0 || y >= nrows())
					continue;
				pr->buf[y][x] = value;
			}
		}else{
			for(int x=x0;x>x1;x--){
				if (x < 0 || x >= ncols())
					continue;
				int y = y0 + dy * (x-x0) / dx;
				if (y < 0 || y >= nrows())
					continue;
				pr->buf[y][x] = value;
			}
		}
	    }else{
		if( y0<y1 ){
			for(int y=y0;y<y1;y++){
				if (y < 0 || y >= nrows())
					continue;
				int x = x0 + dx * (y-y0) / dy;
				if (x < 0 || x >= ncols())
					continue;
				pr->buf[y][x] = value;
			}
		}else{
			for(int y=y0;y>y1;y--){
				if (y < 0 || y >= nrows())
					continue;
				int x = x0 + dx * (y-y0) / dy;
				if (x < 0 || x >= ncols())
					continue;
				pr->buf[y][x] = value;
				}
			}
		}
	    }
	return 1;
	}


/// Treat the picture as a drawing canvas for a
/// time-series graph. There are numvals consecutive
/// points to be plotted over the full width of the
/// picture. The vertical values in yvals are scaled
/// so that minval and maxval are the picture limits.
/// If these are -1, minval and maxval are taken from
/// the data.
	int draw_graph(const int *yvals, const int numvals, 
		const PEL drawshade, int minval = -1, int maxval = -1);

//@}
/// \name Miscellaneous image processing functions
//@{

/// Convolve picture other with a separable point spread function and
/// put the result in this. The horizontal component of the kernel is 
/// given by hsize=m values (h1 h2 .. hm) stored starting at hkernel.
/// The vertical component is given by vsize=n values (v1 v2 .. vn),
/// and the result is to be divided by divisor=k (scaling).
/// All component values must be integers. The kernel is run across the
/// whole picture, and may be preceded by reflection at the boundaries
/// by setting do_symmetric_extension to true. The first argument is
/// not const to allow symmetric extension to be done.
	int sepconv(picture_of_& other, const int hsize, 
		const PEL *hkernel, const int vsize, const PEL *vkernel,
		const bool do_symmetric_extension=false) {
	picture_of_ mid(nrows(),ncols());
	if (do_symmetric_extension)
		other.symmetric_extension(hsize,0);
	mid.hconv(other,hsize,hkernel);
	if (do_symmetric_extension)
		mid.symmetric_extension(0,vsize);
	vconv(mid,vsize,vkernel);
	return 1;
	}

/// Convolve picture other with a horizontal impulse response and put
/// the result in this. The kernel values are stored starting at hkernel
/// and there are hsize of them. You may precede a convolution by a call
/// to symmetric_extension to reflect appropriately at the boundaries if
/// required. You may need to follow a convolution with division by the 
/// sum of the kernel coefficients to normalize the gain.
	PEL hconv(const picture_of_& rb, const int hsize, 
		const PEL *hkernel)
	{
        if ((rowrange.length() != rb.rowrange.length())||
            (colrange.length() != rb.colrange.length()))
                cerr << "warning: picture["<<rowrange.length()<<"]["<<
                        colrange.length()<<
                        "] in a point operation with picture["<<
                        rb.rowrange.length()<<"]["<<
                        rb.colrange.length()<<"]\n";
	PEL accum = 0;
	int outr, inr;
	int outc, inc;
	int hstart = -hsize/2;
	for (outr = rowrange.istart(), inr = rb.rowrange.istart(); 
	    ((outr <= rowrange.iend()) && (inr <= rb.rowrange.iend())); 
	    outr += rowrange.ispacing(), inr += rb.rowrange.ispacing())
	for (outc = colrange.istart(), inc = rb.colrange.istart();
	    ((outc <= colrange.iend()) && (inc <= rb.colrange.iend()));
		outc += colrange.ispacing(), inc += rb.colrange.ispacing()) {
			const PEL *p = hkernel;
			const PEL *q = &rb.pr->buf[inr][inc+hstart];
			PEL pel = 0;
			for (int i = 0; i < hsize; i++)
				pel += *p++ * *q++;
			accum += (pr->buf[outr][outc] = pel);
		}
	return accum;
	}

/// Convolve picture other with a vertical impulse response and put
/// the result in this. The kernel values are stored starting at vkernel
/// and there are vsize of them. You may precede a convolution by a call
/// to symmetric_extension to reflect appropriately at the boundaries if
/// required. You may need to follow a convolution with division by the 
/// sum of the kernel coefficients to normalize the gain.
	PEL vconv(const picture_of_& rb, const int vsize, 
		const PEL *vkernel)
	{
        if ((rowrange.length() != rb.rowrange.length())||
            (colrange.length() != rb.colrange.length()))
                cerr << "warning: picture["<<rowrange.length()<<"]["<<
                        colrange.length()<<
                        "] in a point operation with picture["<<
                        rb.rowrange.length()<<"]["<<
                        rb.colrange.length()<<"]\n";
	PEL accum = 0;
	int outr, inr;
	int outc, inc;
	int vstart = -vsize/2;
	int vinc = (int)(&pr->buf[1][0] - &pr->buf[0][0]);
	for (outr = rowrange.istart(), inr = rb.rowrange.istart(); 
	    ((outr <= rowrange.iend()) && (inr <= rb.rowrange.iend())); 
	    outr += rowrange.ispacing(), inr += rb.rowrange.ispacing())
	for (outc = colrange.istart(), inc = rb.colrange.istart();
	    ((outc <= colrange.iend()) && (inc <= rb.colrange.iend()));
		outc += colrange.ispacing(), inc += rb.colrange.ispacing()) {
			const PEL *p = vkernel;
			const PEL *q = &rb.pr->buf[inr+vstart][inc];
			PEL pel = 0;
			for (int i = 0; i < vsize; i++) {
				pel += *p++ * *q;
				q += vinc;
			}
			accum += (pr->buf[outr][outc] = pel);
		}
	return accum;
	}

/// Convolve picture other with a gaussian filter with standard
/// deviation stddev and put the result in this. The
/// kernel size adapts to the scale of the gaussian. The first argument
/// is not const because gaussfilter does symmetric extension on the
/// source image before convolution.
	int gaussfilter(picture_of_& other, const double stddev) {
	PEL filtvals[511]; // Assume we'll never have longer filter than this
	int filtwidth = (int) (3*stddev*2+1);
	if (!(filtwidth & 1)) filtwidth++;	// Use odd length filters
	// Now make the gaussian
	// cout << "Filter width = " << filtwidth << endl;
	double x,y;
	double variance = stddev*stddev;
	double denom = sqrt(2*3.141592)*stddev;
	int i;
	for (i = 0, x = -3*stddev; i < filtwidth; x += 1.0, i++) {
		y = exp(-x*x/(2.0*variance));
		filtvals[i] = (PEL) (y * 100);
		// The * 100 is so we preserve amplitude resolution for
		// integer types
		}
	// for (i = 0; i < filtwidth; i++)
	// 	cout << filtvals[i] << ' ';
	// cout << endl;
	sepconv(other,filtwidth,filtvals,filtwidth,filtvals,true);
	*this /= (PEL)(100*100*denom*denom);
	return 1;
	}

//@}

/// Waits for given number of milliseconds since last call to wait() on
/// any picture_of_int or colour_picture_of_int. Returns actual number of
/// milliseconds waited (negative if interval has already passed and
/// therefore wait() returned immediately).
	int wait(const int milliseconds);
    };

/// Supports operations on colour pictures via picture_of_int components.

/// colour_picture_of_ maintains four picture_of_int components called mono,
/// red, green and blue.
/// For any colour_picture_of_ either mono alone is used, or red, green	
/// and blue (but not mono) are used. Which of these two cases applies 
/// is determined at construction time and can't be changed.
/// colour_picture_of_ has many functions equivalent
/// to picture_of_'s, except that they
/// operate over all three of red,
/// green and blue if the colour_picture_of_
/// was constructed as colour and over
/// mono if the colour_picture_of_ was
/// constructed as monochrome.
/// Note that colour_picture_of_ omits the
/// following picture_of_int functions: 
/// (1) operator[](const int selector),
///   i.e. You can't do something like
///   a[128][128] = 40;
///   where a is a colour_picture_of_.
///   Instead you must do:
///   if (a.ismonochrome())
///	a.mono[128][128] = 40;
///   else {
///	a.red[128][128] = 40;
///	a.green[128][128] = 40;
///	a.blue[128][128] = 40;
///	}.
/// (2) operator[](const double selector),
///   although map() does work for the
///   the purpose of interpolating one
///   colour_picture_of_ from another.
/// (3) showascolour(const char *comment)
/// and the related reshow functions
///   because a colour_picture_of_ already
///   _is_ colour!
/// (4) mask(const picture_of_& rb,...);
///   if you are masking, you need to
///   do it on monochrome images.
/// (5) operator==(PEL testval)
///   is another operation that only makes
///   sense on monochrome images.
/// (6) bufptr()
///   because colour_picture_of_s have more
///   than one buffer.
/// (7) rowptr(int i)
///   because colour_picture_of_s have more
///   that one row pointer per row.
/// (8) the text constructor,
///   so you must create text via a
///   picture_of_ construction.
/// (9) the matrix functions.
/// (10) The specialized image processing functions for convolutions and the
/// drawing functions (lines, graphs and ramps) are not yet supported.
/// All other functions are supported.
template <typename PEL> class colour_picture_of_ {
public:

    picture_of_<PEL> mono;
    picture_of_<PEL> green;
    picture_of_<PEL> red;	
    picture_of_<PEL> blue;
// Private part of colour_picture_of_int class
private:
    // Following data items are used for sequence processing:
    bool seq;
    int numfiles;
    int currfile;
    bool usingcam;
    int firstcycle;
    char **av;
    colour_picture_of_<PEL> *backup;
    struct clip_internal::CImgDisplay *qxp;
    int puremono;
    colour_picture_of_& assign(const PEL& newval, const char& sel) {
    	if (!puremono) {
    		green.assign(newval,sel);
    		red.assign(newval,sel);
    		blue.assign(newval,sel);
    		}
	else
    		mono.assign(newval,sel);
    	return *this;
    	}


// Assignment for colour_picture with different data type of pels
template <typename OTHERPEL>
	colour_picture_of_& assign(const colour_picture_of_<OTHERPEL>& rb, const char& sel) {
    	if (!puremono) {
		if (rb.ismonochrome()) {
			green.assign(rb.mono,sel);
			red.assign(rb.mono,sel);
			blue.assign(rb.mono,sel);
			}
		else {
    			green.assign(rb.green,sel);
    			red.assign(rb.red,sel);
    			blue.assign(rb.blue,sel);
			}
    		}
	else if (!rb.ismonochrome()) {
		if (sel == '=') {
			mono.assign(rb.green,'=');
			mono.assign(rb.red,'+');
			mono.assign(rb.blue,'+');
			mono.assign(3,'/');		
			}
		else {
			cerr << "Sorry: I can't combine in colour with a colour_picture_of_\n";
			cerr << "that was forced to monochrome during construction.\n";
			cerr << "Make the source monochrome or construct the\n";
			cerr << "destination so that it is colour.\n";
			}
		}
	else
    		mono.assign(rb.mono,sel);
    	return *this;
    	}

// Assignment for colour_picture with same type of pels
	colour_picture_of_& assign(const colour_picture_of_<PEL>& rb) {
    	if (!puremono) {
		if (rb.ismonochrome()) {
			green.assign(rb.mono);
			red.assign(rb.mono);
			blue.assign(rb.mono);
			}
		else {
    			green.assign(rb.green);
    			red.assign(rb.red);
    			blue.assign(rb.blue);
			}
    		}
	else if (!rb.ismonochrome()) {
		mono.average_assign(rb.green,rb.red,rb.blue);
		}
	else
    		mono.assign(rb.mono);
    	return *this;
	}

	PEL sum(const colour_picture_of_& rb, const char& sel) {
	PEL retval;
    	if (puremono != rb.puremono) {
		cerr<< "Sorry: I can't sum up operations between a normal\n";
		cerr<< "colour picture and one that was forced to monochrome\n";
		cerr<< "during construction. Construct the two pictures in\n";
		cerr << "the same way.\n";
		retval = 0;
		}
	else if (!puremono) {
    		retval = green.sum(rb.green,sel);
    		retval += red.sum(rb.red,sel);
    		retval += blue.sum(rb.blue,sel);
		}
	else
    		retval = mono.sum(rb.mono,sel);
    	return retval;
    	}

	// Following constructor sets up dereferenced version
	colour_picture_of_(const irange& range, const colour_picture_of_& source) {
	qxp = 0;	// Never show a dereferenced version
	backup = 0;
	seq = 0;
	puremono = source.puremono;
	mono.dereffrom(range, source.mono);
	red.dereffrom(range, source.red);
	green.dereffrom(range, source.green);
	blue.dereffrom(range, source.blue);
	}


	int construct_or_read_from_file(const char *pnmname,
		const int bord = DEFAULT_BORDER, const int readflag = 1);
// And the part accessible by applications
public:

/// Returns 1 if this picture was constructed as monochrome (e.g.
/// because it was read it from a PGM file), returns 0 if truly colour.
/// USER FUNCTIONS THAT OPERATE AT THE PEL LEVEL ON colour_picture_of_ OBJECTS
/// must always test this. E.g., to read in a picture then put a white
/// dot in the centre: colour_picture_of_ in(argv[1]);
/// int nr = in.nrows(); int nc = in.ncols();
/// if (in.ismonochrome())	in.mono[nr/2][nc/2] = 255;
/// else {	in.green[nr/2][nc/2] = 255;
///	in.red[nr/2][nc/2] = 255;	in.blue[nr/2][nc/2] = 255;
///	}
    int ismonochrome() const 		
	{ return puremono; }

/// Instantiates a colour picture by reading it in to the green, red
/// and blue picture_of_s from a file.  If the file is monochrome, it 
/// is read in to the mono picture_of_ and all subsequent operations are
/// done in mono only. The bord argument works as in picture_of_.
/// CLIP will read the following file formats: PGM, PPM, JPEG, APT and the
/// most common types of BMP.
    colour_picture_of_<PEL>(const char *pnmname,  const int bord=DEFAULT_BORDER) {
	qxp = 0;
	backup = 0;
	seq = 0;
	if (construct_or_read_from_file(pnmname, bord, 0) >= 0)
		return;
	char errmsg[180];
	sprintf(errmsg,"CLIP cannot read\n%s\nas a picture",pnmname);
	clip::picture_of_<PEL> temp(errmsg,0,255);
	int r=temp.nrows();
	int c=temp.ncols();
	puremono = 1;
	mono.initialize(r, c, bord);
	red.fake_initialize(mono);
	green.fake_initialize(mono);
	blue.fake_initialize(mono);
	mono = temp;
	}

	
/// Instantiates a zeroed colour picture of given dimensions. Usually sets up
/// red, green and blue picture_of_s but if forcemono is nonzero, sets
/// up mono only.
    colour_picture_of_(const int r, const int c, const int bord=DEFAULT_BORDER,
		const int forcetomono=0) {
	qxp = 0;
	backup = 0;
	seq = 0;
	if (forcetomono) {
		// cout << "Colour picture being forced to monochrome\n";
		puremono = 1;
		mono.initialize(r, c, bord);
		red.fake_initialize(mono);
		green.fake_initialize(mono);
		blue.fake_initialize(mono);
		return;
		}
	puremono = 0;
	mono.initialize(r, c, bord);
	red.initialize(r, c, bord);
	green.initialize(r, c, bord);
	blue.initialize(r, c, bord);
	}

/// Conventional copy constructor
template <typename OTHERPEL>
    colour_picture_of_ (const colour_picture_of_<OTHERPEL>& source) {
	using clip_internal::pic_rep;
	qxp = 0;
	backup = 0;
	seq = 0;
	puremono = source.puremono;
	// Copy everything except the pic_rep explicitly
	mono.qxp = 0; red.qxp = 0; green.qxp = 0; blue.qxp = 0;
	mono.dereferenced = source.mono.dereferenced;
	mono.rowrange = source.mono.rowrange; 
	mono.colrange = source.mono.colrange;
	red.dereferenced = source.red.dereferenced;
	red.rowrange = source.red.rowrange; 
	red.colrange = source.red.colrange;
	green.dereferenced = source.green.dereferenced;
	green.rowrange = source.green.rowrange; 
	green.colrange = source.green.colrange;
	blue.dereferenced = source.blue.dereferenced;
	blue.rowrange = source.blue.rowrange; 
	blue.colrange = source.blue.colrange;
	// So far haven't copied the pic_reps and don't want to if
	// we are dereferenced
	if (puremono) {
		if (mono.dereferenced) 
			mono.pr = (pic_rep<PEL> *) source.mono.pr;
			// Above cast is ok because never dereference to
			// a different type.
		else 
			mono.pr = new pic_rep<PEL>(*source.mono.pr);
		red.pr = green.pr = blue.pr = mono.pr;
		}
	else {
		if (red.dereferenced) {
			mono.pr = (pic_rep<PEL> *) source.mono.pr;	// Not really necessary
			red.pr = (pic_rep<PEL> *) source.red.pr;
			blue.pr = (pic_rep<PEL> *) source.blue.pr;
			green.pr = (pic_rep<PEL> *) source.green.pr;
			// Above casts are ok because never dereference to
			// a different type.
			}
		else {
			mono.initialize(mono.rowrange.length(),
				mono.colrange.length(),
				source.mono.pr->border);
			red.pr = new pic_rep<PEL>(*source.red.pr);
			green.pr = new pic_rep<PEL>(*source.green.pr);
			blue.pr = new pic_rep<PEL>(*source.blue.pr);
			}
		}
			
	}
// And the special case of constructing from same type that VC doesn't seem to
// deal with automatically.
// template <>
    colour_picture_of_ (const colour_picture_of_<PEL>& source) {
	using clip_internal::pic_rep;
	qxp = 0;
	backup = 0;
	seq = 0;
	puremono = source.puremono;
	// Copy everything except the pic_rep explicitly
	mono.qxp = 0; red.qxp = 0; green.qxp = 0; blue.qxp = 0;
	mono.dereferenced = source.mono.dereferenced;
	mono.rowrange = source.mono.rowrange; 
	mono.colrange = source.mono.colrange;
	red.dereferenced = source.red.dereferenced;
	red.rowrange = source.red.rowrange; 
	red.colrange = source.red.colrange;
	green.dereferenced = source.green.dereferenced;
	green.rowrange = source.green.rowrange; 
	green.colrange = source.green.colrange;
	blue.dereferenced = source.blue.dereferenced;
	blue.rowrange = source.blue.rowrange; 
	blue.colrange = source.blue.colrange;
	// So far haven't copied the pic_reps and don't want to if
	// we are dereferenced
	if (puremono) {
		if (mono.dereferenced) 
			mono.pr = source.mono.pr;
		else 
			mono.pr = new pic_rep<PEL>(*source.mono.pr);
		red.pr = green.pr = blue.pr = mono.pr;
		}
	else {
		if (red.dereferenced) {
			mono.pr = source.mono.pr;	// Not really necessary
			red.pr = source.red.pr;
			blue.pr = source.blue.pr;
			green.pr = source.green.pr;
			}
		else {
			mono.initialize(mono.rowrange.length(),
				mono.colrange.length(),
				source.mono.pr->border);
			red.pr = new pic_rep<PEL>(*source.red.pr);
			green.pr = new pic_rep<PEL>(*source.green.pr);
			blue.pr = new pic_rep<PEL>(*source.blue.pr);
			}
		}
			
	}

/// Creates a colour_picture_of_ from a picture_of_, so loads up mono only.
template <typename OTHERPEL>
    colour_picture_of_ (const picture_of_<OTHERPEL>& source) {	
	qxp = 0;
	backup = 0;
	seq = 0;
	int r, c;
	int bord;
	if (source.dereferenced) {
	    cerr << "Sorry. You cannot construct a colour_picture_of_ from an\n";
	    cerr << "irange-dereferenced picture_of_. Use assignment.\n";
	    exit(1);
	    }
	puremono = 1;
        r = source.rowrange.length();
        c = source.colrange.length();
        bord = source.pr->border;
	mono.initialize(r, c, bord);
	red.fake_initialize(mono);
	green.fake_initialize(mono);
	blue.fake_initialize(mono);
	mono.assign(source,'=');
	}

/// Modify-source constructor that applies callback point function (*f)
/// over source rb. If rb is mono, then so is this picture.
    colour_picture_of_ (const colour_picture_of_<PEL>& source, PEL (*f)(PEL&, const PEL&)) {
	qxp = 0;
	backup = 0;
	seq = 0;
	int forcetomono = source.ismonochrome();
	int r, c;
	int bord;
	if (forcetomono) {
		puremono = 1;
        	r = source.mono.rowrange.length();
        	c = source.mono.colrange.length();
        	bord = source.mono.pr->border;
		mono.initialize(r, c, bord);
		red.fake_initialize(mono);
		green.fake_initialize(mono);
		blue.fake_initialize(mono);
		}
	else {
		puremono = 0;
        	r = source.red.rowrange.length();
        	c = source.red.colrange.length();
        	bord = source.red.pr->border;
		mono.initialize(r, c, bord);
		red.initialize(r, c, bord);
		green.initialize(r, c, bord);
		blue.initialize(r, c, bord);
		}
	point(source, f);
	}

/// Modify-source constructor that applies callback region function (*f)
/// over source rb. If rb is mono, then so is this picture.
    colour_picture_of_(colour_picture_of_<PEL>& source, PEL(*f)(PEL&, const PEL **)) {
 	qxp = 0;
	backup = 0;
	seq = 0;
	int forcetomono = source.ismonochrome();
	int r, c;
	int bord;
	if (forcetomono) {
		puremono = 1;
        	r = source.mono.rowrange.length();
        	c = source.mono.colrange.length();
        	bord = source.mono.pr->border;
		mono.initialize(r, c, bord);
		red.fake_initialize(mono);
		green.fake_initialize(mono);
		blue.fake_initialize(mono);
		}
	else {
		puremono = 0;
        	r = source.red.rowrange.length();
        	c = source.red.colrange.length();
        	bord = source.red.pr->border;
		mono.initialize(r, c, bord);
		red.initialize(r, c, bord);
		green.initialize(r, c, bord);
		blue.initialize(r, c, bord);
		}
	neigh(source, f);
        }

/// Read picture into this object. If this object was constructed as
/// colour and the file is monochrome, will copy the
/// monochrome values to all three of red, green and blue.
/// If this object was constructed as monochrome and the file is colour,
/// read() will refuse to convert.
    int read(const char *ppmname);	

/// Write picture to PPM file (unless it was constructed as monochrome
/// in which case a PGM file will be written.
    int write(const char *ppmname);	

/// Write picture in APT format 0 <= quality <= 100. When quality is
/// 100, coding is lossless
    int aptwrite(const char *aptname, const int quality = 100);

/// Write picture in BMP format. If you must.
    int bmpwrite(const char *bmpname);

/// Destructor
    ~colour_picture_of_() {
	// cerr << "colour picture delete\n";
	if (seq && backup)
		delete backup;
	if (qxp)
		delete qxp;
	if (puremono) {	// We have to set up pseudopics now so that
			// picture_of_ delete works ok
		red.initialize(1,1,0);
		green.initialize(1,1,0);
		blue.initialize(1,1,0);
		}
	}


/// Copies picture_of_ into colour_picture_of_. If this object was
/// constructed as colour, will copy the rb value into red, green and blue.
template <typename OTHERPEL>
    colour_picture_of_<PEL>& operator= (const picture_of_<OTHERPEL>& rb) {
	if (puremono)
		mono.assign(rb,'=');
	else {
		green.assign(rb,'=');
		red.assign(rb,'=');
		blue.assign(rb,'=');
		}
	return *this;
	}




/// Apply callback point function *f() over picture, modifying itself.
    PEL point(PEL (*f)(PEL&)) {
	PEL retval;
	if (puremono)
		return mono.point(f);
	retval = red.point(f);
	retval += green.point(f);
	retval += blue.point(f);
	return retval;
	}

/// Apply callback point function *f() over source picture rb, putting
/// result into *this.
    PEL point(const colour_picture_of_& rb, PEL (*f)(PEL&, const PEL&)) {
	PEL retval;
    	if (!puremono) {
		if (rb.puremono) {
			retval = green.point(rb.mono,f);
			retval += red.point(rb.mono,f);
			retval += blue.point(rb.mono,f);
			}
		else {
    			retval = green.point(rb.green,f);
    			retval += red.point(rb.red,f);
    			retval += blue.point(rb.blue,f);
			}
    		}
	else if (!rb.puremono) {
		cerr << "Sorry: I can't do a colour point operation to a ";
		cerr << "colour_picture_of_\n";
		cerr << "that was forced to monochrome during construction.\n";
		cerr << "Make the source monochrome or construct the\n";
		cerr << "destination so that it is colour.\n";
		retval = 0;
		}
	else
    		retval = mono.point(rb.mono,f);
    	return retval;
	}

/// Apply callback neigh function *f() over source picture rb, putting
/// result into this.
    PEL neigh(colour_picture_of_& rb, PEL(*f)(PEL&, const PEL **)) {
	PEL retval;
    	if (!puremono) {
		if (rb.puremono) {
			retval = green.neigh(rb.mono,f);
			retval += red.neigh(rb.mono,f);
			retval += blue.neigh(rb.mono,f);
			}
		else {
    			retval = green.neigh(rb.green,f);
    			retval += red.neigh(rb.red,f);
    			retval += blue.neigh(rb.blue,f);
			}
    		}
	else if (!rb.puremono) {
		cerr << "Sorry: I can't do a colour neigh operation to a ";
		cerr << "colour_picture_of_\n";
		cerr << "that was forced to monochrome during construction.\n";
		cerr << "Make the source monochrome or construct the\n";
		cerr << "destination so that it is colour.\n";
		retval = 0;
		}
	else
    		retval = mono.neigh(rb.mono,f);
    	return retval;
	}

    
/// \name Assignment operators
//@{
/// Set this picture equal to source picture with different pel data type.
/// This and all assignments
/// work with irange-referenced pictures.
template <typename OTHERPEL>
    colour_picture_of_& operator=(const colour_picture_of_<OTHERPEL>& rb) {
    	return assign(rb, '='); }

/// Set this picture equal to source picture with same pel data type. 
// template <>
    colour_picture_of_& operator=(const colour_picture_of_<PEL>& rb) {
    	return assign(rb); }

/// Set every pel to newval. E.g.irange rg(0,8,256-8);
/// picture_of_ a(256,256); a[rg][rg] = 255;
/// sets every eighth pel horizontally and vertically to white.
    colour_picture_of_& operator=(PEL newval){
    	return assign(newval, '='); }

/// Add source picture to this picture. I.e. each pel value in rb is
/// added to the corresponding pel's value in this.
template <typename OTHERPEL>
    colour_picture_of_& operator+=(const colour_picture_of_<OTHERPEL>& rb) {
    	return assign(rb, '+'); }

/// Add newval to every pel (in all components)
    colour_picture_of_& operator+=(PEL newval){
    	return assign(newval, '+'); }

/// Subtract source picture from this picture. I.e. each pel value
/// in rb is subtracted from the corresponding pel's value in this.
template <typename OTHERPEL>
    colour_picture_of_& operator-=(const colour_picture_of_<OTHERPEL>& rb) {
    	return assign(rb, '-'); }

/// Subtract newval from every pel (in all components)
    colour_picture_of_& operator-=(PEL newval){
    	return assign(newval, '-'); }

/// Multiple source picture with this picture. I.e. each pel value
/// in rb is multiplied with the corresponding pel's value in this.
template <typename OTHERPEL>
    colour_picture_of_& operator*=(const colour_picture_of_<OTHERPEL>& rb) {
    	return assign(rb, '*'); }

/// Multiple every pel by newval (in all components)
    colour_picture_of_& operator*=(PEL newval){
    	return assign(newval, '*'); }


/// If a pel in rb is 0, divides the corresponding pel in this
/// by 1, otherwise divides it by rb pel value.
template <typename OTHERPEL>
    colour_picture_of_& operator/=(const colour_picture_of_<OTHERPEL>& rb) {
    	return assign(rb, '/'); }

/// Divide every pel by newval (unless newval == 0)
    colour_picture_of_& operator/=(PEL newval){
    	return assign(newval, '/'); }

/// The |= operator does a pelwise absolute difference, i.e.
/// destpel = abs(sourcepel-destpel) not a pelwise OR.
template <typename OTHERPEL>
    colour_picture_of_& operator|=(const colour_picture_of_<OTHERPEL>& rb) {
    	return assign(rb, '|'); }

/// Replace every pel with its absolute difference from newval
    colour_picture_of_& operator|=(PEL newval){
    	return assign(newval, '|'); }

/// The ^= operator does a pelwise squared difference, i.e.
/// destpel = sourcepel-destpel; destpel *= destpel;
/// not a pelwise XOR.
template <typename OTHERPEL>
    colour_picture_of_& operator^=(const colour_picture_of_<OTHERPEL>& rb) {
    	return assign(rb, '^'); }

/// Replace every pel with its squared difference from newval
    colour_picture_of_& operator^=(PEL newval){
    	return assign(newval, '^'); }

/// The &= operator does a pelwise conditional replacement, i.e.
/// if (sourcepel < 0) destpel=destpel; else destpel = sourcepel;
template <typename OTHERPEL>
    colour_picture_of_& operator&=(const colour_picture_of_<OTHERPEL>& rb) {
    	return assign(rb, '&'); }

// Following two operators were included in earlier versions of
// clip for completeness, but they really have no value!
    colour_picture_of_& operator&=(PEL newval){
    	return assign(newval, '&'); }

//@}

/// \name Operators combining pictures to give PEL results
//@{

/// When a picture is added to a picture the result is a PEL which is the
/// sum of all pelwise additions.
    PEL operator+(const colour_picture_of_& rb) {
    	return sum(rb, '+'); }

/// Return the sum of all pelwise subtractions.
    PEL operator-(const colour_picture_of_& rb) {
    	return sum(rb, '-'); }

/// Return the sum of all pelwise multiplications. Can convolve a picture with
/// an operator by making the operator into a picture, then multiply operator
/// and a moving block of the picture together.
    PEL operator*(const colour_picture_of_& rb) {
    	return sum(rb, '*'); }

/// Return the sum of all pelwise divisions. Uses divisor of 1 for any zero
/// pel in rb.
    PEL operator/(const colour_picture_of_& rb) {
    	return sum(rb, '/'); }

/// Return the sum of all pelwise absolute differences. Can be used for block matching.
    PEL operator|(const colour_picture_of_& rb) {
    	return sum(rb, '|'); }

/// Return the sum of all pelwise squared differences. Can be used for block matching.
    PEL operator^(const colour_picture_of_& rb) {
    	return sum(rb, '^'); }

//@}

/// \name Dereferencing operators
/// The following operator[]s are used to dereference (i.e. index) pictures.
/// They must always be appended to a colour_picture_of_ in pairs, with the
/// first [] specifying a row value or range and the second [] specifying
/// a column value or range. 

//@{

/// When dereferenced (i.e. indexed) with irange operators, a subimage or
/// block of the picture is returned as a picture. E.g. irange rg(0,1,7);
/// colour_picture_of_ pic(256,256); pic[rg][rg] = 255; sets the top left 8x8 square 
/// in pic to white.
    colour_picture_of_ operator[](const irange& selector) {
	colour_picture_of_ m(selector,*this);
	return m;
	}

/// We often want to paste one picture on another, centred at a particular
/// location. Can do this with irange derefencings, but this way is easier.
    int paste_centred_on (const picture_of_<PEL>& from,
		const int x, const int y) {
	irange r(0,1,from.nrows()-1);
	irange c(0,1,from.ncols()-1);
	r -= from.nrows()/2;
	c -= from.ncols()/2;
	r += y;
	c += x;
	if (puremono)
		mono[r][c].assign(from,'=');
	else {
		red[r][c].assign(from,'=');
		green[r][c].assign(from,'=');
		blue[r][c].assign(from,'=');
		}
	// operator[](r).operator[](c).assign(from,'=');
	return 0;
	}

/// As above, for colour_picture_of_ pasting
    int paste_centred_on (const colour_picture_of_<PEL>& from,
		const int x, const int y) {
	irange r(0,1,from.nrows()-1);
	irange c(0,1,from.ncols()-1);
	r -= from.nrows()/2;
	c -= from.ncols()/2;
	r += y;
	c += x;
	operator[](r).operator[](c).assign(from,'=');
	return 0;
	}

//@}

/// \name Picture display functions
/// Use these functions to show or reshow the current picture on the screen

//@{

/// Display the picture on the screen in its own window, with
/// the comment string as title. If normalize is true values are scaled
/// to [0,255] for display. If false, it is assumed they already are in this range.
    int show(const char *comment, const bool normalize = 1);

/// Display picture without a specific title
    int show(const bool n = 1) { return show("CLIP image", n); }

/// Display without decorations at Nx,Ny. If parent==0, the coords are
/// screen coordinates, and the window receives events (e.g. pointx())
/// as normal. If parent is a displayed picture, the coords are relative to
/// that picture and a child window is displayed that does not receive events.
    int showat(const int Nx, const int Ny, const bool normalize,
			const picture_of_<PEL> *parent);

/// Display without decorations, with colour_picture_of_ as parent
    int showat(const int Nx, const int Ny, const bool normalize,
			const colour_picture_of_<PEL> *parent);

/// Display without decorations at screen coords Nx, Ny
	int showat(const int Nx, const int Ny, const bool normalize) {
		return showat(Nx,Ny,normalize,(const picture_of_<PEL> *)0); }

/// Display without decorations at screen coords Nx, Ny
 	int showat(const int Nx, const int Ny) {
		return showat(Nx,Ny,1,(const picture_of_<PEL> *)0); }

/// Redisplay the picture (presumably because it has changed) and change
/// title to the comment string.
    int reshow(const char *comment, const bool normalize = 1);

/// Redisplay the picture without altering the title.
    int reshow(const bool normalize = 1);

/// Redisplay without normalization. Equivalent to reshow(comment,0);
    int fastreshow(const char *comment);

/// Redisplay without normalization. Equivalent to reshow(0);
    int fastreshow();

/// Change title of the displayed pic.
    int rename(const char *comment);

/// Stop showing the picture
    int unshow();

/// Returns non-zero if the user has left-clicked in the close box of the
/// window where the picture is being displayed.
    int closerequested();		

/// Return current pointer (mouse) x coord relative to left of picture window    
	int pointx();

/// Return current pointer (mouse) y coord relative to top of picture window    
    int pointy();

/// Return non-zero if the mouse's left button is pressed within the picture's window.
    int pointleftbutton();
	
/// Returns non-zero if mid button down
	int pointmiddlebutton();
	
/// Returns non-zero if right button down
    int pointrightbutton();	
	
/// Returns non-zero if the shift key is held down while a mouse button is down.
    int pointshiftbutton();		

/// Returns keycode of any key held down    
	int key();

/// Clears keystroke so that it will not be reprocessed.
    int clearkey();	

/// Inspect a shown picture. Until the right button is clicked inside
/// the picture, the program pauses and allows the user to
/// left-click within the picture and inspect the values at that point.	
    int inspect() {
	if (!qxp) return -1;
	char infostring[64];
	while (pointrightbutton());	// From any previous rightbutton press
	int prevx = -1;
	int prevy = -1;
	while(1) {
		if (pointrightbutton()) {
			// No longer do a rename back to what was before
			return 1;
			}
		if (!pointleftbutton())
			continue;
		int x = pointx();
		int y = pointy();
		if ((x == prevx)&&(y == prevy))
			continue;
		prevx = x; prevy = y;
		if ((x < 0)||(x >= ncols())||(y < 0)||(y >= nrows())) {
			// Beyond picture limits
			rename("pointer outside pic");
			continue;		// Back to start of loop
			}
		// Else pointer is within picture limits
		if (puremono) {
			float pelvalue = (float) mono.pr->buf[y][x];
			sprintf(infostring,"value at (%d,%d) = %f",
				x,y,pelvalue);
			}
		else {
			float redvalue = (float) red.pr->buf[y][x];
			float greenvalue = (float) green.pr->buf[y][x];
			float bluevalue = (float) blue.pr->buf[y][x];
			sprintf(infostring,"value at (%d,%d) = %f,%f,%f",
				x,y,redvalue,greenvalue,bluevalue);
			}
		rename(infostring);	// Update header
		}
	}


//@}

/// Set values in the border of the picture by copying the values on
/// the edges of the picture outwards.
    void dc_pad() {
	if (puremono)
		mono.dc_pad();
	else {
		red.dc_pad();
		green.dc_pad();
		blue.dc_pad();
		}
	}


/// Set values in the border of the picture to zero.
    void zero_border() { 
	if (puremono)
		mono.zero_border();
	else {
		red.zero_border();
		green.zero_border();
		blue.zero_border();
		}
	}

/// Return number of rows in picture				
    int nrows() const { return mono.nrows(); }

/// Return number of columns in picture
    int ncols() const { return mono.ncols(); }

/// Return width of border
    int bordersize() { return mono.bordersize(); }

/// Return maximum value in picture
    PEL max() const {
    	if (puremono)
    		return mono.max();
    	PEL r = red.max(); 
	PEL g = green.max(); 
	PEL b = blue.max();
    	if (g > r) r = g;
    	if (b > r) r = b;
    	return r;
    	}

/// Clip picture at given maximum value
    void max(const PEL clipval) {
	if (puremono)
		mono.max(clipval);
	else {
		green.max(clipval);
		red.max(clipval);
		blue.max(clipval);
		}
	}

/// Return minimum value in picture
PEL min() const { 
    	if (puremono)
    		return mono.min();
    	PEL r = red.min(); 
	PEL g = green.min(); 
	PEL b = blue.min();
    	if (g < r) r = g;
    	if (b < r) r = b;
    	return r;
    	}

/// Clip picture at given minimum value
    void min(const PEL clipval) {
	if (puremono)
		mono.min(clipval);
	else {
		green.min(clipval);
		red.min(clipval);
		blue.min(clipval);
		}
	}

 
/// Return total of all pel values
PEL total() const {
	if (puremono)
		return mono.total();
	return (red.total() + green.total() + blue.total());
	}


/// Reranges or normalizes picture values to fall in given range
/// If scale_even_if_inside_range is false (default), then if all
/// values are current between the give newlow and newhigh, they
/// are unaltered. Otherwise, values are scaled so that the new
/// max() will be newhigh and the new min() will be newmin.
    int rerange(const PEL newlow, const PEL newhigh, 
    			const bool scale_even_if_inside_range = 0) {
    	if (newlow == newhigh) {
		*this = newhigh;
		return 0;
		}
    	PEL mx = max();
	PEL mn = min();
	if ((mx == newhigh) && (mn == newlow))
		return 0;
	if (!scale_even_if_inside_range && (mx <= newhigh) && (mn >= newlow))
		return 0;
	if (mx == mn) {
		if (mx > 0)	*this = newhigh;
		else		*this = newlow;
		return 0;
		}
	PEL newrange = newhigh - newlow;
	PEL oldrange = mx - mn;
	if (puremono) {
	    for (int i = 0; i < nrows(); i++) {
		PEL *pel = mono.pr->buf[i];
		for (int j = 0; j < ncols(); j++) {
			*pel = newlow + (newrange*(*pel-mn))/oldrange;
			pel++;
			}
		}
	    }
	else {
	    for (int i = 0; i < nrows(); i++) {
		PEL *r = red.pr->buf[i];
		PEL *g = green.pr->buf[i];
		PEL *b = blue.pr->buf[i];
		for (int j = 0; j < ncols(); j++) {
			*r = newlow + (newrange*(*r-mn))/oldrange;
			*g = newlow + (newrange*(*g-mn))/oldrange;
			*b = newlow + (newrange*(*b-mn))/oldrange;
			r++, g++, b++;
			}
		}
	    }
	return 1;
	}

/// \name Sequence processing facilities
/// Sequence processing facilities for "standard" CLIP programs.
/// Many CLIP programs consist of a processing loop that cycles over
/// input pictures from a sequence of named files or from a camera.

//@{

/// Sequence processing constructor that interprets main()'s parameters. 
/// If the command line had no arguments, tries to attach a camera. 
/// If that fails, generates a default picture. 
/// If the command line had arguments, interprets these as picture
/// files and sets up for reading these in sequence.
/// See also member function next()
colour_picture_of_<PEL>(int argc, char *argv[], 
		const int def_nrows=240, const int def_ncols=320, 
		const int bord=DEFAULT_BORDER) {
	seq = 1;	// Flag to say sequence processing
	numfiles = argc - 1;
	currfile = 1;
	usingcam = 0;
	firstcycle = 1;
	av = argv;
	qxp = 0;
	backup = 0;
	if (!numfiles) {
		// If there are no arguments, try to open a camera and, if
		// that fails, create a default image
		puremono = 0;
		mono.initialize(def_nrows, def_ncols, bord);
		red.initialize(def_nrows, def_ncols, bord);
		green.initialize(def_nrows, def_ncols, bord);
		blue.initialize(def_nrows, def_ncols, bord);
		if (attach_camera() < 0) {
			red.vramp();
			blue.hramp();
			green = 128;
			clip::picture_of_<PEL> text("Default image",0,128,3);
			green.paste_centred_on(text,
						def_ncols/2,def_nrows/2);
			}
		else {
			request_frame();
			usingcam = 1;
			}
		return;
		}
	// Otherwise, open the first file
	if (construct_or_read_from_file(argv[1], bord, 0) >= 0) {
		seq = 1;
		if (numfiles == 1)
			// Make and save a backup for fast reread
			backup = new colour_picture_of_(*this);
		return;
		}
	qxp = 0;
	char errmsg[180];
	sprintf(errmsg,"CLIP cannot read\n%s\nas a picture",argv[1]);
	clip::picture_of_<PEL> temp(errmsg,0,255);
	int r=temp.nrows();
	int c=temp.ncols();
	puremono = 1;
	mono.initialize(r, c, bord);
	red.fake_initialize(mono);
	green.fake_initialize(mono);
	blue.fake_initialize(mono);
	mono = temp;
	if (numfiles == 1)
		// Make and save a backup for fast reread
		backup = new colour_picture_of_(*this);
	}

/// Loads the next picture in the sequence into this picture that was
/// constructed with the sequence processing constructor. Returns 1 if
/// the picture has been seen before (if cycling back to the start of
/// the argument list) and 0 if it has not. Returns -1 if the picture
/// was not constructed using the sequence processing constructor.
int next()
	{
	if (!seq)	// Not set up to do sequence processing
		{
		cerr << "next() called for a picture not constructed for ";
		cerr << "sequence processing." << endl;
		return -1;
		}
	if (usingcam) {
		while(!frame_delivered())
			;
		request_frame();
		return 1;	// Not seen this picture before
		}
	if (!numfiles){		// Using default picture
		// Remake it, just in case it was overwritten
		red.vramp();
		blue.hramp();
		green = 128;
		clip::picture_of_<PEL> text("Default image",0,128,3);
		green.paste_centred_on(text,ncols()/2,nrows()/2);
		return 0;
		}
	if (numfiles == 1){	// We made a backup of this one to save
				// rereading every cycle
		if (puremono)
			mono = backup->mono;
		else {
			red = backup->red;
			blue = backup->blue;
			green = backup->green;
			}
		return 0;
		}
	if (currfile > numfiles) {
		currfile = 1;	// Back to start
		firstcycle = 0;
		}
	read(av[currfile++]);
	return firstcycle;
	}

//@}
	
		
/// The map operators provide a fast way of mapping a bilinearly interpolated
/// value from one picture into another. transparency is a number from 0 (overwrite) 
/// to 255 (preserve almost all of what was there originally). So with two
/// colour_picture_of_s's source and dest, dest.map(50,76,source,5.3,99.99,0);
/// is exactly equivalent to dest.comp[76][50] = source.comp[99.99][5.3]; for
/// comp = {red, green, blue} or comp = {mono}, but much faster.
    void map(const int myx,const int myy,const colour_picture_of_& from,
	const double fromx,const double fromy, const int transparency);


/// This map operator is even faster than the double version above. This
/// time fromx and fromy are both 256 times the actual values desired, so
/// dest.map(50,76,source,64,384,0); is the same as 
/// dest.map(50,76,source,0.25,1.5,0); but faster. Rarely used.
    void map(const int myx,const int myy,const colour_picture_of_& from,
	const int fromx,const int fromy, const int transparency);

/// \name Functions for interfacing to a local camera
//@{

/// If a camera exists, attach it to this picture. Returns negative value 
/// if the camera doesn't exist, could not be initialized, is already 
/// attached to another picture or the current picture is too small to 
/// have a camera attached.
	int attach_camera();

/// Initiate the request of a frame from the camera attached to this 
/// picture. Returns immediately; you must check frame_delivered() for 
/// arrival of the frame.
	int request_frame();

/// This version of request_frame() only works with Philips cameras on 
/// Linux. Will be extended to other cameras and platforms when possible.
	int request_frame(int shutterspeed, int gain);

/// Immediately returns 0 if a requested frame is not yet available. If 
/// the frame is available does a buffer switch and returns 1. Because of 
/// the buffer switch within the function, immediately before the call the 
/// picture has its old contents; immediately after it has its the frame.
	int frame_delivered();

/// Call to detach the camera from the picture. You may need to do this 
/// in order to attach the camera to a different picture.
	int detach_camera();
//@}

/// Waits for given number of milliseconds since last call to wait() on
/// any picture_of_ or colour_picture_of_. Returns actual number of
/// milliseconds waited (negative if interval has already passed and
/// therefore wait() returned immediately).
	int wait(const int milliseconds);

// Have to deal with friends differently because of MSC VC++ problem
#if ( defined(_MSC_VER)) 
	friend class picture_of_<unsigned char>;
	friend class picture_of_<char>;
	friend class picture_of_<unsigned short>;
	friend class picture_of_<short>;
	friend class picture_of_<unsigned int>;
	friend class picture_of_<int>;
	friend class picture_of_<unsigned long>;
	friend class picture_of_<long>;
	friend class picture_of_<float>;
	friend class picture_of_<double>;
	friend class picture_of_<long double>;
	friend class colour_picture_of_<unsigned char>;
	friend class colour_picture_of_<char>;
	friend class colour_picture_of_<unsigned short>;
	friend class colour_picture_of_<short>;
	friend class colour_picture_of_<unsigned int>;
	friend class colour_picture_of_<int>;
	friend class colour_picture_of_<unsigned long>;
	friend class colour_picture_of_<long>;
	friend class colour_picture_of_<float>;
	friend class colour_picture_of_<double>;
	friend class colour_picture_of_<long double>;
#else
	template <class APEL> friend class picture_of_;
	template <class APEL> friend class colour_picture_of_;
#endif
    };

// Finally the stuff on double_dereferenced pictures.
// Not to be used by applications.
namespace clip_internal {
template <typename PEL> class double_dereferenced_picture_of_ {
	PEL *abovep;
	PEL *belowp;
	double a;
public:
	double_dereferenced_picture_of_<PEL>(const double& r, picture_of_<PEL> *source) {
		int top = (int) r;
		if (r < 0)
			top--;
		a = r - top;
		abovep = source->pr->buf[top++]; 
		belowp = source->pr->buf[top]; 
		}
	PEL operator[](const double& c) {
		int left = (int) c;
		if (c < 0)
			left--;
		double b = c - left;
		int y00 = abovep[left];
		int y10 = belowp[left++];
		int y01 = abovep[left];
		int y11 = belowp[left];
		return (PEL)(a*b*(y11 - y01 - y10 + y00) + a*(y10 - y00) + b*(y01 - y00) + y00 + 0.5);
		}
	};
}	// namespace clip_internal

/// Supports parameter control through interaction on picture windows

/// parampic is a particular kind of colour_picture_of_ providing easy facilities
/// for a programmer to display one- and two-dimensional sliders via CLIP.
/// It can be supplied with a background picture, or use default ramps, then 
/// it displays a hairline slider over the background. The range of the 
/// slider is set during construction, and thereafter, a user clicking in 
/// the picture moves the slider to the appropriate scaled value. The program 
/// can then easily read off the current value of the parameter(s) 
/// represented by the parampic. It is assumed that the parameter read 
/// functions are called often when input is expected because sensing of pointer
/// (mouse) state and updating are only done when the functions are called.
template <typename PEL> class parampic_of_{
	colour_picture_of_<PEL> *bg;
	colour_picture_of_<PEL> *shown;
	bool bg_is_mine;
	int lastx, lasty;
	int type;
	double scalex, scaley;
	double lowx, lowy;
	double x, y;	// Current values
	void paintvalues(double xval, double yval) {
		char text[128];
		if (type == 3)
			sprintf(text,"(%.3f,%.3f)",xval,yval);
		else if (type == 2)
			sprintf(text,"(%.3f)",yval);
		else
			sprintf(text,"(%.3f)",xval);
		picture_of_<PEL> t(text,255,-1);
		irange r(0,1,t.nrows()-1);
		irange c(0,1,t.ncols()-1);
		if (shown->ismonochrome())
			shown->mono[r][c] &= t;
		else {
			shown->red[r][c] &= t;
			shown->green[r][c] &= t;
			shown->blue[r][c] &= t;
			}
		}
	void paintcross(int x, int y) {
		irange yr(0,2,shown->nrows()-1);
		irange xr(0,2,shown->ncols()-1);
		irange thickness(0,0,0);
		if (shown->ismonochrome()) {
			if (type & 1){
				shown->mono[yr][thickness+x] = 0;
				shown->mono[yr+1][thickness+x] = 255;
				}
			if (type & 2){
				shown->mono[thickness+y][xr] = 0;
				shown->mono[thickness+y][xr+1] = 255;
				}
		}
		else {
			if (type & 1) {
				shown->red[yr][thickness+x] = 0;
				shown->green[yr][thickness+x] = 0;
				shown->blue[yr][thickness+x] = 0;
				shown->red[yr+1][thickness+x] = 255;
				shown->green[yr+1][thickness+x] = 255;
				shown->blue[yr+1][thickness+x] = 255;
			}
			if (type & 2) {
				shown->red[thickness+y][xr] = 0;
				shown->green[thickness+y][xr] = 0;
				shown->blue[thickness+y][xr] = 0;
				shown->red[thickness+y][xr+1] = 255;
				shown->green[thickness+y][xr+1] = 255;
				shown->blue[thickness+y][xr+1] = 255;
			}
		}
	}
	// init is called by the constructors
	void init (colour_picture_of_<PEL> *background, 
		const double xstart, const double xend,
		const double xfirst, const double ystart, 
		const double yend, const double yfirst,
		const char *name) {
		if (ystart == yend)
			type = 1;	// Horizonal slider only
		else if (xstart == xend)
			type = 2;	// Vertical slider only
		else
			type = 3;	// 2D
		if (!background) {
			if (type == 1) {
				bg = new colour_picture_of_<PEL>(32,256,0,1);
				bg->mono.hramp();
				}
			else if (type == 2) {
				bg = new colour_picture_of_<PEL>(256,64,0,1);
				// Vertical sliders are wider than horiz
				// ones are tall so there's room for
				// a caption
				bg->mono.vramp();
				}
			else {
				bg = new colour_picture_of_<PEL>(256,256,0,0);
				bg->red.hramp();	
				bg->green.vramp();
				bg->blue = 128;
				}
			bg_is_mine = true;
			}
		else {
			bg = background;
			bg_is_mine = false;
			}
		shown = new colour_picture_of_<PEL>(*bg);
		lowx = xstart;
		lowy = ystart;
		double xrange = xend-xstart;
		double yrange = yend-ystart;
		scalex = xrange/bg->ncols();
		scaley = yrange/bg->nrows();
		x = xfirst;
		y = yfirst;
		lastx = (int)((xfirst-xstart)/scalex);
		lasty = (int)((yfirst-ystart)/scaley);
		paintcross(lastx,lasty);
		paintvalues(x,y);
		shown->show(name);
	}
public:

/// Create a parameter picture with the given background or a default if the
/// background pointer argument = 0, with X,Y control ranges as specified.
/// The X parameter value can vary between xstart and xend, with initial value
/// xfirst; the Y parameter value can vary between ystart and yend, with initial
/// value yfirst. If xstart == xend, the slider is 1D and vertical;
/// if ystart == yend, the slider is 1D and horizontal. Otherwise, both X and Y
/// can be varied.
	parampic_of_(colour_picture_of_<PEL> *background, 
			const double xstart, const double xend,
			const double xfirst, const double ystart, 
			const double yend, const double yfirst,
			const char *name) {
		init(background,xstart,xend,xfirst,ystart,yend,yfirst,name);
		}

/// Create the most common kind of parameter picture. 1D in X with default
/// background. 
	parampic_of_(const double xstart, const double xend, const double xfirst,
			const char *name) {
		init(0,xstart,xend,xfirst,0,0,0,name);
		}

/// Destructor
	~parampic_of_() {
		delete shown;
		if (bg_is_mine)
			delete bg;
	}

/// Get current X and Y values into xout and yout. Returns 1 if there has
/// been a change, 0 otherwise.
	int XY(double &xout, double &yout) {
		xout = x; yout = y;
		if (!shown->pointleftbutton())
			return 0;
		if ((lastx == shown->pointx())&&(lasty == shown->pointy()))
			return 0;
		lastx = shown->pointx();
		lasty = shown->pointy();
		xout = x = ((double)lastx)*scalex + lowx;
		yout = y = ((double)lasty)*scaley + lowy;
		*shown = *bg;
		paintcross(lastx,lasty);
		paintvalues(x,y);
		shown->reshow();
		return 1;
	}

/// Get current X value into xout. Returns 1 if there has
/// been a change, 0 otherwise.
	int X(double &xout) {
		xout = x;
		if (!shown->pointleftbutton())
			return 0;
		if (lastx == shown->pointx())
			return 0;
		lastx = shown->pointx();
		xout = x = ((double)lastx)*scalex + lowx;
		*shown = *bg;
		paintcross(lastx,lasty);
		paintvalues(x,y);
		shown->reshow();
		return 1;
	}
/// Get current Y value into yout. Returns 1 if there has
/// been a change, 0 otherwise.
	int Y(double &yout) {
		yout = y;
		if (!shown->pointleftbutton())
			return 0;
		if (lasty == shown->pointy())
			return 0;
		lasty = shown->pointy();
		yout = y = ((double)lasty)*scaley + lowy;
		*shown = *bg;
		paintcross(lastx,lasty);
		paintvalues(x,y);
		shown->reshow();
		return 1;
	}
};

/// \name Input and output of pictures as text

//@{

/// Output monochrome picture as newline delimited rows of space 
/// delimited cols.
template<class T>
std::ostream& operator<<(std::ostream& os, picture_of_<T>& pic)
{
	for(int i = 0; i < pic.nrows(); ++i) {
		for(int j = 0; j < pic.ncols(); ++j)
			os << pic[i][j] << ' ';
		os << std::endl;
		}
	return os;
}

/// Input monochrome picture from text
template<class T>
std::istream& operator>>(std::istream& is, picture_of_<T>& pic)
{
	for(int i = 0; i < pic.nrows(); ++i) {
		for(int j = 0; j < pic.ncols(); ++j)
			is >> pic[i][j];
		}
	return is;
}


/// Output colour picture as double-newline delimited planes of
/// newline delimited rows of space delimited cols.
template<class T>
std::ostream& operator<<(std::ostream& os, colour_picture_of_<T>& pic)
{
	for(int i = 0; i < pic.nrows(); ++i) {
		for(int j = 0; j < pic.ncols(); ++j)
			os << pic.red[i][j] << ' ';
		os << std::endl;
		}
	os << std::endl;
	for(int i = 0; i < pic.nrows(); ++i) {
		for(int j = 0; j < pic.ncols(); ++j)
			os << pic.green[i][j] << ' ';
		os << std::endl;
		}
	os << std::endl;
	for(int i = 0; i < pic.nrows(); ++i) {
		for(int j = 0; j < pic.ncols(); ++j)
			os << pic.blue[i][j] << ' ';
		os << std::endl;
		}
	return os;
}

/// Input colour picture from text.
template<class T>
std::istream& operator>>(std::istream& is, colour_picture_of_<T>& pic)
{
	for(int i = 0; i < pic.nrows(); ++i) {
		for(int j = 0; j < pic.ncols(); ++j)
			is >> pic.red[i][j];
		}
	for(int i = 0; i < pic.nrows(); ++i) {
		for(int j = 0; j < pic.ncols(); ++j)
			is >> pic.green[i][j];
		}
	for(int i = 0; i < pic.nrows(); ++i) {
		for(int j = 0; j < pic.ncols(); ++j)
			is >> pic.blue[i][j];
		}
	return is;
}

//@}


/* CLIP_E */

//
// New platform-independent display facilities via part of the CImg
// library. See the copyright notice for the library below.
// Modifications to this part of CImg.h are scattered throughout.
// John Robinson. September 2004.
//
/*------------------------------------------------------------------------------------------------------
  
  File        : CImg.h
  Description : The C++ Template Image Processing Library
  Author      : David Tschumperl�
  Institution : Initially developped at ODYSSEE Lab, INRIA Sophia Antipolis.
  
  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.
  
  ----------------------------------------------------------------------------------------------------*/

/*
  Detect and set CImg configuration flags

  
  If compilation flags are not adapted to your system, you can
  set them to correct values, before including the file CImg.h.

  
  -------------------------------------------------------------*/

} // namespace clip
// Overcome VisualC++ 6.0 and DMC compilers namespace bug
#if ( defined(_MSC_VER) && _MSC_VER<=1200 ) || defined(__DMC__)
#define std
#endif

// Autodetection of the current OS
#ifndef cimg_OS
#if defined(NO_OS_SUPPORT)
// Compiler flag to omit operating system support
#define cimg_OS           -1
#ifndef cimg_display_type
#define cimg_display_type  0
#endif
#elif defined(sun) || defined(__sun)        
// Solaris
#define cimg_OS            0
#ifndef cimg_display_type
#define cimg_display_type  1
#endif
#elif defined(linux) || defined(__linux) || defined(__CYGWIN__)
// Linux
#define cimg_OS            1
#ifndef cimg_display_type
#define cimg_display_type  1
#endif
#ifndef cimg_color_terminal
#define cimg_color_terminal
#endif
#elif defined(_WIN32) || defined(__WIN32__)
// Windows
#define cimg_OS            2
#ifndef cimg_display_type
#define cimg_display_type  2
#endif
#elif defined(__MACOSX__) || defined(__APPLE__)
// Mac OS X
#define cimg_OS            3
#ifndef cimg_display_type
#define cimg_display_type  1
#endif
#else 
// Other configurations
#define cimg_OS           -1
#ifndef cimg_display_type
#define cimg_display_type  0
#endif
#endif
#endif
// Debug flag. Set to 0 to remove dynamic debug messages.
#ifndef cimg_debug
#define cimg_debug         1
#endif

// Architecture-dependent includes
#if cimg_OS==-1
// Don't include anything
#elif cimg_OS!=2
#include <sys/time.h>
#include <unistd.h>
#else
#ifndef _WIN32_WINNT
#define _WIN32_WINNT 0x0400     // Overcome bug in macro definition in winuser.h (allow the use of TrackMouseEvent() )
#endif
#include <windows.h>
#endif
#if cimg_display_type==1
#include <X11/Xlib.h>
#include <X11/Xutil.h>
#include <pthread.h>
namespace clip {
namespace clip_internal {
static pthread_mutex_t* cimg_mutex = NULL;
static pthread_t*       cimg_event_thread = NULL;
static CImgDisplay*     cimg_wins[1024];
static Display*         cimg_displayX11 = NULL;
static volatile unsigned int    cimg_nb_wins = 0;
static volatile bool    cimg_thread_finished = false;
static unsigned int             cimg_bits = 0;
static GC*              cimg_gc = NULL;
}	// namespace clip_internal
} 	// namespace clip
#endif

#ifdef min
#undef min
#undef max
#undef abs
#endif



#if cimg_OS==0 || cimg_OS==1
  // Keycodes for Unix & Solaris
  const unsigned int keyESC        = 9;
  const unsigned int keyF1         = 67;
  const unsigned int keyF2         = 68;
  const unsigned int keyF3         = 69;
  const unsigned int keyF4         = 70;
  const unsigned int keyF5         = 71;
  const unsigned int keyF6         = 72;
  const unsigned int keyF7         = 73;
  const unsigned int keyF8         = 74;
  const unsigned int keyF9         = 75;
  const unsigned int keyF10        = 76;
  const unsigned int keyF11        = 95;
  const unsigned int keyF12        = 96;
  const unsigned int keyPAUSE      = 110;
  const unsigned int key1          = 10;
  const unsigned int key2          = 11;
  const unsigned int key3          = 12;
  const unsigned int key4          = 13;
  const unsigned int key5          = 14;
  const unsigned int key6          = 15;
  const unsigned int key7          = 16;
  const unsigned int key8          = 17;
  const unsigned int key9          = 18;
  const unsigned int key0          = 19;
  const unsigned int keyBACKSPACE  = 22;
  const unsigned int keyINSERT     = 106;
  const unsigned int keyHOME       = 97;
  const unsigned int keyPAGEUP     = 99;
  const unsigned int keyTAB        = 23;
  const unsigned int keyQ          = 24;
  const unsigned int keyW          = 25;
  const unsigned int keyE          = 26;
  const unsigned int keyR          = 27;
  const unsigned int keyT          = 28;
  const unsigned int keyY          = 29;
  const unsigned int keyU          = 30;
  const unsigned int keyI          = 31;
  const unsigned int keyO          = 32;
  const unsigned int keyP          = 33;
  const unsigned int keyDELETE     = 107;
  const unsigned int keyEND        = 103;
  const unsigned int keyPAGEDOWN   = 105;
  const unsigned int keyCAPSLOCK   = 66;
  const unsigned int keyA          = 38;
  const unsigned int keyS          = 39;
  const unsigned int keyD          = 40;
  const unsigned int keyF          = 41;
  const unsigned int keyG          = 42;
  const unsigned int keyH          = 43;
  const unsigned int keyJ          = 44;
  const unsigned int keyK          = 45;
  const unsigned int keyL          = 46;
  const unsigned int keyENTER      = 36;
  const unsigned int keySHIFTLEFT  = 50;
  const unsigned int keyZ          = 52;
  const unsigned int keyX          = 53;
  const unsigned int keyC          = 54;
  const unsigned int keyV          = 55;
  const unsigned int keyB          = 56;
  const unsigned int keyN          = 57;
  const unsigned int keyM          = 58;
  const unsigned int keySHIFTRIGHT = 62;
  const unsigned int keyARROWUP    = 98;
  const unsigned int keyCTRLLEFT   = 37;
  const unsigned int keyAPPLEFT    = 115;
  const unsigned int keySPACE      = 65;
  const unsigned int keyALTGR      = 113;
  const unsigned int keyAPPRIGHT   = 116;
  const unsigned int keyMENU       = 117;
  const unsigned int keyCTRLRIGHT  = 109;
  const unsigned int keyARROWLEFT  = 100;
  const unsigned int keyARROWDOWN  = 104;
  const unsigned int keyARROWRIGHT = 102;

#elif cimg_OS==3 
  // Keycodes for MacOS X
/*
Well the following codes are those in the CImg package, but
my Mac generates a different set, so those follow this comment
  const unsigned int keyESC        = 61;
  const unsigned int keyF1         = 130;
  const unsigned int keyF2         = 128;
  const unsigned int keyF3         = 107;
  const unsigned int keyF4         = 126;
  const unsigned int keyF5         = 104;
  const unsigned int keyF6         = 105;
  const unsigned int keyF7         = 106;
  const unsigned int keyF8         = 108;
  const unsigned int keyF9         = 109;
  const unsigned int keyF10        = 117;
  const unsigned int keyF11        = 111;
  const unsigned int keyF12        = 119;
  const unsigned int keyPAUSE      = 110;
  const unsigned int key1          = 26;
  const unsigned int key2          = 27;
  const unsigned int key3          = 28;
  const unsigned int key4          = 29;
  const unsigned int key5          = 31;
  const unsigned int key6          = 30;
  const unsigned int key7          = 34;
  const unsigned int key8          = 36;
  const unsigned int key9          = 33;
  const unsigned int key0          = 37;
  const unsigned int keyBACKSPACE  = 59;
  const unsigned int keyINSERT     = 122;
  const unsigned int keyHOME       = 123;
  const unsigned int keyPAGEUP     = 124;
  const u