[MINC-users] coordinate question

[Steve M. Robbins]

On Mon, Jul 10, 2006 at 04:16:22PM -0400, Andrew Janke wrote:

> 2. transform your world co-ordinate using something like octave.
> 
> 
> I don't know of a MINC specific command that will transform an input
> point from the C/L. The other option is to convert your world
> co-ordinate to a tagfile and then transform the tag file.

Here's a simple program I wrote years ago that transforms points.
Do with it what you like.

-Steve


/* Apply an MNI general transform to one or more points.
 */

// TODO: add facility to invert transform

#include <iostream>
#include <string>


extern "C" {
#   include <volume_io.h>
}
#include <ParseArgv.h>

using namespace std;




/* Read one point per line from in, transform, and
 * write on out.
 */
void transform_points( std::istream& in,
		       General_transform& transform,
		       std::ostream& out )
{
    while( !in.eof() ) {
	Real x,y,z;

	in >> x;
	in >> y;
	in >> z;

	if ( in.eof() )  break;

	Real tx,ty,tz;
	general_transform_point( &transform,
				 x, y, z,
				 &tx, &ty, &tz );

	out << tx << " " << ty << " " << tz << endl;
    }
}
	

const char* transform_type_name[] = {
    "linear", "thin plate spline", "user-defined",
    "concatenated", "grid"
};


/* Diagnostics
 */
void display_transform_info( std::ostream& os,
			     General_transform* t,
			     std::string prefix = "" )
{
    os << prefix << "Type of transform: ";
    if ( t->type > GRID_TRANSFORM ) {
	os << "type " << t->type << " [UNKNOWN]";
    } else {
	os << transform_type_name[t->type];
    }

    os << endl
       << prefix << "Invert flag: " << t->inverse_flag << endl;

    if ( t->type == CONCATENATED_TRANSFORM ) {
	int n = get_n_concated_transforms(t);
	for( int i = 0; i < n; ++i ) {
	    display_transform_info( os, get_nth_general_transform(t,i),
				    prefix + "    " );
	}
    }
}



int main( int ac, char* av[] ) 
{
    static int verbose = 0;

    ArgvInfo argTable[] = {
	{ "-verbose", 
	  ARGV_CONSTANT, (char*)1, (char*)&verbose,
	  "emit diagnostic info" },
	{ NULL, ARGV_END, NULL, NULL, NULL }
    };

    
    if ( ParseArgv( &ac, av, argTable, 0 ) || ac != 2 ) {
	cerr << "usage: " << av[0] 
	     << " transformation"
	     << endl;
	return 1;
    }

    General_transform transform;
    if ( input_transform_file( av[1], &transform ) != OK ) {
	cerr << av[0] << ": cannot read transform file: " 
	     << av[1] << endl;
	return 2;
    }

    if (verbose) 
	display_transform_info( std::cerr, &transform );

    transform_points( std::cin, transform, std::cout );

    return 0;
}