squirmgrid.cpp

本程序模拟细胞的自我繁殖

			a=b;
		}
		// up the right
		for(i=length-1;i>=0;i--)
		{
			b = CreateNewCell(start_x+(width+1)*RADIUS*2.0F,start_y+(i+1)*RADIUS*2.0F,'a',S);
			b->BondTo(a);
			a=b;
		}
		// along the top
		for(i=width;i>=0;i--)
		{
			b = CreateNewCell(start_x+(i+1)*RADIUS*2.0F,start_y,'a',S);
			b->BondTo(a);
			a=b;
		}
		// link to the start of the loop
		a->BondTo(start);
	}

    // create some random cells around the place
    for(int i=0;i<n;i++)
        CreateRandomCell();

    this->initial_population = n;
}

float sqr(float x) { return x*x; }

void SquirmGrid::DoTimeStep()
{
    RecomputeVelocitiesAndReact();
    MoveCells();
    EdgeEffect();
    if( this->do_flood && (this->iterations%this->FLOOD_PERIOD)==0 && this->iterations>0 )
        DoFlood();
    this->iterations++;
}

SquirmCell* SquirmGrid::CreateRandomCell()
{
    // pick a random position
    C2dVector loc;
    loc.x = RandomX();
    loc.y = RandomY();
	return CreateNewCell(loc,"abcdef"[RandomN(6)],0);//RandomN(4)!=0?'c':(RandomN(2)==0?'a':'b'),0);
}

SquirmCell* SquirmGrid::CreateRandomCellOfType(char type)
{
    // pick a random  position
    C2dVector loc;
    loc.x = RandomX();
    loc.y = RandomY();
    return CreateNewCell(loc.x,loc.y,type,0);
}

SquirmCell* SquirmGrid::CreateNewCell(C2dVector loc,char type,int state)
{
    // only proceed if this location is empty
    if(Valid(loc))
    {
        SquirmCell *cell = new SquirmCell(loc,type,state);

        this->PutCell(cell);
        this->cell_list.AddTail(cell);

        return cell;
    }
    else
    {
        // this position is taken!
        SquirmError("CreateNewCell : slot taken!");
        return NULL;
    }
}

void SquirmGrid::PutCell(SquirmCell *cell)
{
	// put the cell into a slot
	int cx = (int)floor(cell->location.x/slot_size);
	int cy = (int)floor(cell->location.y/slot_size);
    cx = max(0,cx);
    cx = min(cx,this->slots_X);
    cy = max(0,cy);
    cy = min(cy,this->slots_Y);
	cell_grid[cx][cy].addOccupant(cell);
}

void SquirmGrid::RemoveCell(SquirmCell *cell)
{
	// clear the slot it is in
	int cx = (int)floor(cell->location.x/slot_size);
	int cy = (int)floor(cell->location.y/slot_size);
    cx = max(0,cx);
    cx = min(cx,this->slots_X);
    cy = max(0,cy);
    cy = min(cy,this->slots_Y);
	cell_grid[cx][cy].removeOccupant(cell);
}
void SquirmGrid::Draw(CDC* pDC,float scale)
{
    //if(this->iterations%10!=0) return;

    if(!this->created) return;

    // outline the area
    pDC->Rectangle(0,0,(int)(this->size.x*scale)+OFFSET*2,
		(int)(this->size.y*scale)+OFFSET*2);

    if(this->probe)
    {
        // draw the probe
        pDC->FillSolidRect((int)(probe->location.x*scale)-2+OFFSET,
		    (int)(probe->location.y*scale)-2+OFFSET,6,6,
		    probe->GetColour());
    }

    // draw the other atoms
    POSITION pos = this->cell_list.GetHeadPosition();
    SquirmCell *cell,*other;
    bool first=true;
    while(pos)
    {
        cell = this->cell_list.GetNext(pos);
        // draw the cell as a filled box or a pixel if window is small
        if(scale>=0.5F)
            pDC->FillSolidRect((int)(cell->location.x*scale)-1+OFFSET,
			  (int)(cell->location.y*scale)-1+OFFSET,3,3,
			  cell->GetColour());
        else
            pDC->SetPixelV((int)(cell->location.x*scale)+OFFSET,
			  (int)(cell->location.y*scale)+OFFSET,cell->GetColour());
        // optionally, draw as a RADIUS circle
        if(true && cell->GetState()!=0)
        {
            CRect r(int(cell->location.x*scale-scale*RADIUS)+OFFSET,int(cell->location.y*scale-scale*RADIUS)+OFFSET,
                int(cell->location.x*scale+scale*RADIUS)+OFFSET,int(cell->location.y*scale+scale*RADIUS)+OFFSET);
            CPen pen(PS_SOLID,1,cell->GetColour()),*old_pen;
            CBrush brush,*old_brush;
            brush.CreateSolidBrush(cell->GetColour());
            old_brush = pDC->SelectObject(&brush);
            old_pen = pDC->SelectObject(&pen);
            pDC->Ellipse(&r);
            pDC->SelectObject(old_pen);
            pDC->SelectObject(old_brush);
        }
		// draw the cell's state
		if(false)
		{
			CString text;
			text.Format("%d",cell->GetState());
			pDC->TextOut((int)(cell->location.x*scale)+OFFSET,
				(int)(cell->location.y*scale)+OFFSET,text);
		}
        // draw the cell's bonds as lines
		if(true) {
			POSITION bonds_pos = cell->GetBondedCells().GetHeadPosition();
			while(bonds_pos)
			{
				other = cell->GetBondedCells().GetNext(bonds_pos);
				pDC->MoveTo(int(cell->location.x*scale)+OFFSET,
					int(cell->location.y*scale)+OFFSET);
				pDC->LineTo(int(other->location.x*scale)+OFFSET,
					int(other->location.y*scale)+OFFSET);
			}
		}
    }
}

void SquirmGrid::DoFlood()
{
    // for all the atoms in one part of the area:
    // - remove all their bonds
    // - set their state to 0
    SquirmCell *cell;

    float left[4]={0,this->GetSize().x/2.0F,this->GetSize().x/2.0F,0};
    float top[4]={0,0,this->GetSize().y/2.0F,this->GetSize().y/2.0F};


    POSITION pos = this->cell_list.GetHeadPosition();
	while(pos)
	{
		cell = this->cell_list.GetNext(pos);
        if( cell->location.x>=left[this->flood_sector] && 
            cell->location.x<=left[this->flood_sector]+this->GetSize().x/2.0F &&
            cell->location.y>=top[this->flood_sector] && 
            cell->location.y<=top[this->flood_sector]+this->GetSize().y/2.0F )
        {
			// break all of its bonds
			cell->BreakAllBonds();
			// set its state to 0
			cell->SetState(0);
            // assign a random velocity
            cell->AssignRandomVelocity();
        }
    }

    // do the flood on the other side next time
    this->flood_sector = (this->flood_sector+1)%4;
}


float SquirmGrid::RandomX()
{
    return rand()*size.x/RAND_MAX;
}

float SquirmGrid::RandomY()
{
    return rand()*size.y/RAND_MAX;
}

int SquirmGrid::RandomN(int n)
{
    return (rand()*100+rand()) % n;
}

/*CString SquirmGrid::GetReport(CPoint p,int scale)
{
    int x,y;
    x = p.x/scale;
    y = p.y/scale;
    if(!Valid(x,y) || cell_grid[x][y].queryEmpty())
        return "";
    else
    {
        CString report;
        report.Format("%c%d",cell_grid[x][y].getOccupant()->GetType(),
            cell_grid[x][y].getOccupant()->GetState());
        return report;
    } 
}*/

void SquirmGrid::GetAllWithinRadius(C2dVector loc,float r,
											  SquirmCellList &dest)
{
	dest.RemoveAll();
	// which slot is at the centre?
	int cx = (int)floor(loc.x/slot_size);
	int cy = (int)floor(loc.y/slot_size);
    float r2 = r*r;
	if(cx>=0 && cy>=0 && cx<this->slots_X && cy<this->slots_Y)
	{
		int search_slots = (int)ceil(r / slot_size);
		int x,y;
        SquirmCellSlot *slot;
		SquirmCell *cell;
        POSITION pos;
		for(x=max(0,cx-search_slots);x<=min(this->slots_X-1,cx+search_slots);x++)
		{
			for(y=max(0,cy-search_slots);y<=min(this->slots_Y-1,cy+search_slots);y++)
			{
                slot = &cell_grid[x][y];
				if(!slot->queryEmpty())
				{
					// now do a more expensive test
					pos = slot->occupants.GetHeadPosition();
					while(pos)
					{
						cell = slot->occupants.GetNext(pos);
                        if(C2dVector::Dist2(cell->location,loc)<r2)
							dest.AddTail(cell);
					}
				}
			}
		}
	}
	/* case below if required
    else {
		AfxMessageBox("Out of area in GetAllWithinRadius!");
		ASSERT(false);
		return;
	}*/
}

void SquirmGrid::InitSearchArea(int cx,int cy,int search_slots,SquirmCellList*& central_cells,
                                SquirmCellList &nearby_cells)
{
    // assuming the cx,cy is valid
    central_cells = &this->cell_grid[cx][cy].occupants;
    nearby_cells.RemoveAll();
    int i,j;
    for(i=max(0,cx-search_slots);i<=min(this->slots_X-1,cx+search_slots);i++)
        for(j=max(0,cy-search_slots);j<=min(this->slots_Y-1,cy+search_slots);j++)
            nearby_cells.AddTail(&this->cell_grid[i][j].occupants);
}

void SquirmGrid::RecomputeVelocitiesAndReact()
{
    const float PHYS_RANGE = RADIUS*2.0F; // how far out do the physical forces extend?
    const float PHYS_RANGE2 = PHYS_RANGE*PHYS_RANGE;

    const float REACTION_RANGE = RADIUS*2.5F; // how far out can reactions extend?
    const float REACTION_RANGE2 = REACTION_RANGE*REACTION_RANGE;

	const float MAX_RANGE = max(PHYS_RANGE2,REACTION_RANGE2); // no effect acts beyond this distance

    // find the search radius distance in terms of the number of slots
	int search_slots = (int) ceil( MAX_RANGE / this->slot_size);

    // for each slot in the grid, get the list of cells within the neighbourhood
    // use this list to recompute the velocities of the cells in the central slot

    int cx,cy; // central square 
    SquirmCellList *central_cells; // the cells in the central square
    SquirmCellList nearby_cells; // the cells in the search area (including central_cells)