How to Write New Image Functions

Overview

Imfit comes with a variety of 2D surface-brightness functions (“image functions”). But it is designed to make it relatively simple to add new image functions, to describe objects or substructures which are not well modeled by the existing functions. This is done by writing additional C++ code and re-compiling the various Imfit programs (imfit, imfit-mcmc, makeimage) to include the new function or functions.

This page provides some guidelines on how to write new image functions (see also the “Rolling Your Own Image Functions” section of the Imfit manual, from which this page is excerpted). Even if you’re not very familiar with C++, it should hopefully not be too difficult to write a new image function, since it can be done by copying and modifying one of the existing image functions.

A new image function is written in C++ as a subclass of the FunctionObject base class, which is declared and defined in the function_object.h and function_object.cpp source-code files (in the function_objects/ subdirectory of the source-code distribution).

A new image-function class should provide its own implementation of (i.e., “override”) the following public FunctionObject methods, which are defined as virtual methods in the base class:

  • The class constructor – in most cases, the code for this can be copied from any of the existing FunctionObject subclasses, unless some special extra initialization is needed.

  • Setup() – this is used by the calling program to supply the current set of function parameters (including the (x0,y0) pixel values for the center of the function set) prior to determining intensity values for individual pixels. The parameter values in the input array should be assigned to the appropriate data members of the class. This is also a convenient place to do any general calculations which depend on the parameter values but don’t depend on the exact pixel (x,y) values.

  • GetValue() – this is used by the calling program to obtain the surface brightness for a given pixel location (x,y). In existing FunctionObject subclasses, this method often calls other (private) methods to handle details of the calculation.

  • GetClassShortName() – this is a class function which returns the short version of the class name as a string.

The new class should also redefine the following internal class constants:

  • N_PARAMS — the number of input parameters (excluding the central pixel coordinates);

  • PARAM_LABELS — a vector of string labels for the input parameters;

  • FUNCTION_NAME — a short string describing the function;

  • className — a string (no spaces allowed) giving the official name of the function.

The add_functions.cpp file should then be updated by:

  • including the header file for the new class;

  • adding 2 lines to the PopulateFactoryMap() function to add the ability to create an instance of the new class. (Look for the comment line “// ADD CODE FOR NEW FUNCTIONS HERE” in the file.)

Finally, the name of the C++ implementation file for the new class should be added to the SConstruct file to ensure it gets included in the compilation; the easiest thing is to add the file’s name (without the .cpp suffix) to the multi-line functionobject_obj_string string definition. (Or look for the comment line “# ADD CODE FOR NEW FUNCTIONS HERE” in the SConstruct file.)

Existing examples of FunctionObject subclasses can be found in the function_objects/ subdirectory of the source-code distribution, and are the best place to look in order to get a better sense of how to implement new FunctionObject subclasses.

A Simple Example

To demonstrate the basics of writing a new image function, we’ll modify the existing Gaussian class to make a class called NewMoffat, which produces an elliptical structure with a Moffat radial profile. (Such a function already exists in Imfit under the name Moffat, so this is really a redundant exercise.)

Three basic changes to the existing func_gauss.h/cpp files are needed:

  1. Change the class name (in this case, from “Gaussian” to “NewMoffat”);

  2. Change the code which actually computes the function;

  3. Add, rename, and delete class data members to accommodate the new algorithm.

1. Create and Edit the Header File

Assuming we’re in the Imfit source code directory, cd to the function_objects/ subdirectory, copy the header file func_gaussian.h to func_new-moffat.h, and edit the file to change the following lines:

#define CLASS_SHORT_NAME "Gaussian"

class Gaussian : public FunctionObject

Gaussian( );

to

#define CLASS_SHORT_NAME "NewMoffat"

class NewMoffat : public FunctionObject

NewMoffat( );

2. Create and Edit the Class File

Copy the file func_gaussian.cpp to func_new-moffat.cpp.

  1. Change the following lines in the beginning of the file:

#include "func_gaussian.h"

const int N_PARAMS = 4;
const char PARAM_LABELS[][20] = {"PA", "ell", "I_0", "sigma"};
const char  PARAM_UNITS[][30] = {"deg (CCW from +y axis)", "", "counts/pixel", "pixels"};
const char FUNCTION_NAME[] = "Gaussian function";

to

#include "func_new-moffat.h"

const int N_PARAMS = 5;
const char PARAM_LABELS[][20] = {"PA", "ell", "I_0", "fwhm", "beta"};
const char  PARAM_UNITS[][30] = {"deg (CCW from +y axis)", "", "counts/pixel",
                            "pixels", ""};
const char FUNCTION_NAME[] = "Moffat function";

B. In the remainder of the file, change all references to the class name from Gaussian to NewMoffat (e.g., Gaussian::Setup becomes NewMoffat::Setup).

C. Change the Setup method. Here, you’ll need to change how the input parameter array is converted into individual parameters, and do any useful pre-computations (i.e., computations that depend on the parameter values, but not on individual pixel values or values derived from the latter, like radius).

Change

PA = params[0 + offsetIndex];
ell = params[1 + offsetIndex];
I_0 = params[2 + offsetIndex];
sigma = params[3 + offsetIndex];

to

PA = params[0 + offsetIndex];
ell = params[1 + offsetIndex];
I_0 = params[2 + offsetIndex];
fwhm = params[3 + offsetIndex];
beta = params[4 + offsetIndex];

Then, at the end of the method, replaced this line

twosigma_squared = 2.0 * sigma*sigma;

with this (which computes the “alpha” parameter of the Moffat function)

double exponent = pow(2.0, 1.0/beta);
alpha = 0.5*fwhm/sqrt(exponent = 1.0);

  1. Changes to the CalculateIntensity method:

Although it is the public method GetValue which is called by other parts of the program, we don’t actually need to change the current version of that method in this example. The code in the original Gaussian version of GetValue converts pixel positions to a scaled radius value, given input values for the center, ellipticity, and position angle, and then calls the private method CalculateIntensity to determine the intensity as a function of the radius. Since we’re still assuming a perfectly elliptical shape, we can keep the existing code. (GetValue also includes possible pixel subsampling, which is useful for cases where intensity changes rapidly one scales of a single pixel; we’ll apply a simple modification for the Moffat function later on.)

So in this case we actually implement the details of the new function’s algorithm in CalculateIntensity. Replace the original version of that method with the following:

double NewMoffat::CalculateIntensity( double r )
{
  double  scaledR, denominator;

  scaledR = r / alpha;
  denominator = pow((1.0 + scaledR*scaledR), beta);
  return (I_0 / denominator);
}

  1. Changes to the CalculateSubsamples method:

Although pixel subsampling is performed in the GetValues method, the determination of whether or not to actually *do** the subsampling – and how much of it to do – is determined in CalcualteSubsamples.

For the Gaussian function, subsampling can be useful happen when r < 1 and sigma < 1. The equivalent for the Moffat function would be r < 1 and alpha < 1, so change the line in CalculateSubsamples that says

if ((sigma <= 1.0) && (r <= 1.0))

to say

if ((alpha <= 1.0) && (r <= 1.0))

At this point, most of the work is done. We only need to update the code in add_functions.cpp so it knows about the new function and update the SConstruct file so that the new function is included in the compilation.

Adding Parameter Unit Definitions

Imfit version 1.9 added the possibility of including strings describing the intended units for certain parameters of the image functions. This is purely cosmetic – the image functions ignore them – but it can be useful for the user as a reminder of what individual parameters mean. (Example: the Exponential image function has the units label “deg (CCW from +y axis)” for the PA (position angle) parameter, “counts/pixel” for I_0, and “pixels” for the scale length parameter h.)

You do not have to provide this information in any image functions you write. If you wish to, you need to properly define the PARAM_UNITS constant at the top of the source file (see 2.A above for an example); assign the unit-description strings (the empty string – “” – is used for an parameters without units) to the parameterUnits vector in the constructor (see any of the standard image functions for an example of how to do this), and set the parameterUnitsExist data member to true (again, see the constructor of any of the standard image functions for examples).

Other Potential Issues

If your new image function has an analytic expression for the total flux, then you might consider overriding the CanCalculateTotalFlux method to return true and then overriding the TotalFlux method so that it calculates and returns the total flux. (The default is to let makeimage estimate the total flux numerically, by generating a large image using the image function and summing all the pixel values.)

If your new image function is meant to represent the image background (as in the case of the built-in functions FlatSky and TiltedSkyPlane), then you may not want makeimage trying to calculate the “total flux” for the component. In this case, you should add a line to your class’s constructor which sets the data member isBackground to true (as is done in, e.g., func_flatsky.cpp).